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Abstract 

This  thesis  develops  and  evaluates  a  number  of  new  concepts  and  tools  for  the  analysis  of 
signals  using  Malvar  wavelets  (lapped  orthogonal  transforms).  Windowing,  often  employed  as  a 
spectral  estimation  technique,  can  result  in  irreparable  distortions  in  the  transformed  signal.  By 
utilizing  the  Malvar  wavelet,  any  signal  distortion  resulting  from  the  transformation  can  be  elimi¬ 
nated  or  cancelled  during  reconstruction.  This  is  accomplished  by  placing  conditions  on  the  window 
and  the  basis  function  and  then  incorporating  the  window  into  the  orthonormal  representation. 

Paradigms  for  both  a  complex-valued  and  a  real-valued  Malvar  wavelet  are  summarized.  The 
algorithms  for  the  wavelets  were  implemented  in  the  DOD  standard  language,  Ada.  The  code  was 
verified  to  ensure  perfect  reconstruction  could  be  obtained  and  experiments  were  performed  using 
the  wavelet  programs.  Various  compression  techniques  were  investigated  and  evaluated  using  the 
Malvar  wavelet  in  both  homomorphic  and  non-homomorphic  systems.  The  Malvar  wavelet  has  the 
unique  ability  to  overlap  adjacent  windows  without  increasing  the  number  of  transform  coefficients. 
Various  amounts  of  window  overlap  were  investigated  to  determine  if  there  is  an  optimal  amount 
which  should  be  used.  In  addition,  the  real- valued  basis  function  was  used  to  attempt  real- valued 
deconvolution.  It  was  found  that  the  read-valued  Malvar  wavelet,  with  just  one  point  of  overlap 
for  each  window,  provided  better  or  equally  as  good  reconstruction  of  signals  as  most  of  the  more 
complex  systems.  This  same  read-valued  baisis  function  could  be  used  to  perform  deconvolution,  if 
the  original  signal  has  certain  symmetries. 
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FILTERING,  CODING,  AND  COMPRESSION 
WITH  MALVAR  WAVELETS 


I.  Introduction 


Background 

The  United  States  has  a  great  need  for  improved  digital  signal  processing  techniques  in  order 
to  reduce  the  burden  on  overloaded  communications  satellites  and  to  reap  more  information  from 
satellite  imagery. 

Data  Compression:  One  of  the  largest  problems  during  Operation  Desert  Storm  was  the  lack 
of  communication  channels  available  to  field  and  command  units.  Virtually  100%  of  the  military 
ultra-high  frequency  and  super-high  frequency  satellite  communication  channels  were  consumed 
by  logistics,  administrative,  and  intelligence  data  during  the  conflict  (5:82).  Over  90%  of  the 
communications  into  and  out  of  the  theater  went  over  communication  satellites,  however  over  20% 
of  the  information  had  to  be  sent  via  commercial  satellites  (36:4). 

Compressing  speech  requires  a  three  way  trade-off  among  the  goals  of  preserving  intelligi¬ 
bility  and  quality,  limiting  the  bit  rate,  and  minimizing  computational  complexity  (23:225).  All 
three  of  these  objectives  must  be  considered  when  chocsing  the  best  compression  algorithm.  Trans¬ 
form  coding  is  one  of  the  major  low  bit  rate  speech  coding  techniques  being  investigated  by  the 
military  (38). 

Imagery:  Space  imagery  often  requires  adaptive  restoration  to  deblur  out-of-focus  imagery 
because  it  is  generally  impossible  to  merely  retake  a  given  photograph  under  better  conditions  (12). 
“Digital  image  restoration  is  useful  whenever  needed  information  in  an  image  is  hidden  by  some 
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type  of  degradation.  The  classic  statement  of  the  image  restoration  problem  is:  given  a  noisy  and 
blurred  image,  find  an  estimate  of  the  ideal  image  using  a  priori  information  about  the  blur,  noise, 
and  the  ideal  image”  (12).  The  research  for  this  thesis  deals  only  with  one  dimensional  signals, 
however  many  of  the  techniques  can  easily  be  adapted  to  two  dimensional  images.  The  Malvar 
wavelet  is  becoming  increasingly  important  to  all  fields  of  digital  signal  processing  and  h»s  been 
successfully  used  with  an  imaging  technique  called  motion  estimation  (41). 

Whenever  a  signal  is  transmitted  through  any  medium  it  becomes  convolved  with  a  function 
related  to  that  medium.  Speech  is  a  signal  that  is  a  convolution  of  the  speaker’s  voice  pitch  and 
the  sound  of  the  word  made  by  the  vocal  tract.  A  satellite  image  can  be  represented  by  the  original 
image  convolved  with  its  optical  transfer  function  (or  blurring  function,  which  is  a  mathematical 
representation  of  the  atmosphere).  A  seismic  recording  can  be  described  as  the  original  signal 
convolved  with  impulses,  which  appear  as  echoes  in  the  seismic  event.  The  term  deconvolution 
indicates  the  technique  used  to  remove  the  transfer  function  from  a  given  signal.  This  includes  the 
process  of  debluring  images  (12)  or  restoring  old  acoustic  recordings  (27). 

Processing  Technique:  The  Lapped  Orthogonal  Transform  (LOT)  is  an  important  signal 
processing  tool  that  has  been  used  to  filter  and  code  both  speech  (26)  and  image  (17)  (41)  signals. 
The  LOT  was  first  introduced  by  H.S.  Malvar  and  D.H.  Staelin,  in  1988,  in  order  to  reduce  errors  in 
image  reconstruction  (17).  The  complete  derivation  and  implementation  of  the  LOT  is  outlined  in 
Chapter  2.  The  LOT  has  also  been  termed  the  Malvar  wavelet,  by  Yves  Meyer  (19),  because  of  the 
algorithm’s  importance  to  wavelet  analysis;  complimenting  wavelet  packet  analysis.  Ville  proposed 
two  types  of  signal  analysis  in  1947.  He  explained  that  “we  can  either:  first  cut  a  signal  into  slices 
(in  time)  with  a  switch;  then  pass  these  different  slices  through  a  system  of  filters  to  analyze  them. 
Or  we  can:  first  filter  different  frequency  bands;  then  cut  these  bands  into  slices  (in  time)  to  study 
their  energy  variations  (19).”  Yves  Meyer  has  explained  that  the  first  approach  leads  us  to  “Malvar 
wavelets”  and  the  second  approach  leads  to  “wavelet  packets”  (19).  These  two  algorithms  are  the 
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underlying  foundation  for  wavelet  analysis  (19)  (29).  The  terms  “Lapped  Orthogonal  Transform” 
and  “Malvar  wavelet”  will  be  used  interchangeably  throughout  this  paper. 

Problem  Statement 

Using  the  Malvar  wavelet,  investigate  new  ways  in  which  this  signal  processing  technique  may 
be  used  to  filter,  code,  and  compress  digital  signals. 

Research  Goals 

The  research  for  this  thesis  effort  was  essentially  a  set  of  experiments  all  related  to  the  Malvar 
wavelet  (LOT)  analysis.  Several  goals  were  established  which  guided  the  research  for  this  thesis. 
The  research  goals  are  as  follows: 

1.  Investigate  the  use  of  a  complex- valued  basis  with  the  Malvar  wavelets.  Can  orthogonality  be 
preserved  and  can  the  algorithm  be  efficiently  coded?  Can  this  wavelet  be  successfully  used 
to  perform  data  compression? 

2.  Investigate  the  use  of  the  Lapped  Orthogonal  Transform  in  a  homomorphic  filter  for  speech 
processing.  Determine  if  there  is  a  best  method  for  data  compression. 

3.  Investigate  the  use  of  overlap  with  real-valued  Malvar  wavelets.  Is  there  a  correlation  between 
the  amount  of  overlap  and  the  amount  of  compression? 

4.  Investigate  the  use  of  the  sinusoidal  basis  function  (as  described  in  Chapter  II  and  used  with 
the  Sinusoidal  Malvar  wavelet)  in  a  homomorphic  filter  to  perform  signal  deconvolution.  Can 
real-valued  basis  functions  be  useful  for  signal  deconvolution? 

Outline 

Chapter  II  expands  on  the  background  and  theory  used  for  this  research  effort.  It  includes 
development  of  homomorphic  filtering  and  the  lapped  transform.  A  complex- valued  basis  function 


1-3 


has  many  inherent  advantages  and  Chapter  III  reports  the  effort  to  use  an  exponential  basis  func¬ 
tion  with  the  Malvar  wavelet.  Chapter  IV  presents  the  results  of  the  lapped  transform  evaluated 
with  numerous  compression  techniques,  including  both  homomorphic  and  non-homomorphic  filters. 
Chapter  V  is  concerned  with  finding  the  optimal  amount  of  overlap  within  the  real-valued  Mal¬ 
var  wavelet.  Homomorphic  deconvolution  was  investigated  and  the  results  of  this  work  constitute 
Chapter  VI.  Finally,  a  summary  of  key  results,  and  a  recommendation  for  follow-on  research  areas, 
are  outlined  in  Chapter  VII. 
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II.  Background 


Introduction 

Two  important  signal  processing  techniques  must  be  well-defined  for  this  research  effort: 

1.  The  Homomorphic  Filter 

2.  The  Malvar  Wavelet  (Lapped  Orthogonal  Transform) 

This  chapter  will  review  homomorphic  systems,  the  cepstrum  (a  specific  case  of  the  homo¬ 
morphic  filter),  transforms,  their  measurements  of  performance,  and  finally  the  lapped  transform. 

Homomorphic  Filters 

Homomorphic  filters  are  non-linear  digital  signal  processing  (DSP)  systems  that  are  very 
important  to  speech  and  image  processing  (22).  Homomorphic  systems  were  originally  described 
by  Dr.  Alan  Oppenheim  (21).  Oppenheim  showed  that,  in  its  most  general  form,  a  homomorphic 
system  can  be  represented  by  algebraically  linear  (homomorphic)  mappings  between  the  input 
and  output  signals  (22:771).  This  system  can  be  broken  down  into  three  separate  homomorphic 
system  (as  shown  in  figure  2.1);  a  characteristic  system  (D),  a  linear  system  (L),  and  the  inverse 
characteristic  system  (D(inv))  (22:772). 


Figure  2.1.  Canonic  representation  of  homomorphic  systems 

The  cepstrum  is  a  popular  dsp  tool  that  is  a  special  case  of  a  homomorphic  system.  The 
cepstrum  implements  the  Fourier  transform  as  its  characteristic  system  and  the  natural  logarithm 
as  its  linear  system  (27). 
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Cepstrum 

Bogert,  Healy,  and  Tukey  first  defined  the  cepstrum  by  paraphrasing  the  term  “spectrum”, 
because  cepstrum  techniques  perform  a  further  spectral  analysis  on  a  frequency  spectrum  (4).  Using 
the  same  naming  convention,  terms  such  as  “quefrency”  and  “littering”  correspond  to  frequency 
and  filtering  in  the  cepstral  domain  (6:1429).  Two  versions  of  the  cepstrum  are  used  in  signal 
processing  applications;  the  complex  cepstrum  and  the  power  cepstrum. 

The  power  cepstrum  is  defined  as:  “the  power  spectrum  of  the  logarithm  of  the  power  spec¬ 
trum”  (27).  This  is  denoted  as: 

c,  =  |P[log(|X(/)|2)]|2  (2.1) 

where  |X(/)|2  is  the  power  spectrum  of  the  signal  and  T  is  the  Fourier  transform.  Note  that  the 
final  squaring  operation  is  really  unnecessary  and  is  often  not  used  (27).  The  independent  variable 
of  the  cepstrum  is  termed  “quefrency”  (4),  which  is  a  time  signal  in  terms  of  frequency  (units 
of  time  which  are  treated  like  frequency).  The  Fourier  transform  of  a  spectrum  yields  time  on 
the  independent  axis.  A  large  time  value  (high  quefrency)  represents  a  rapid  fluctuation  in  the 
spectrum  and  a  low  quefrency  represents  a  slow  fluctuation  (27).  Thus  the  cepstrum  can  space 
apart  these  fluctuations  in  the  the  frequency  of  the  signal  and  ultimately  separate  two  convolved 
signals  in  this  quefrency  domain. 

Speech  processing  has  made  great  use  of  the  power  cepstrum’s  ability  to  separate  signals. 
Voiced  speech  is  made  up  of  two  distinct  signals;  the  voice  pitch  (p(x))  and  the  sound  of  the  word 
from  the  vocal  tract  <\{x)).  These  two  signals  are  convolved  together  to  make  coherent  speech; 
s(x)  =  p(x)  *  v(x)  (23).  Because  of  the  convolution  property  of  the  Fourier  transform,  the  power 
spectrum  of  the  signal  changes  the  cor  solution  into  multiplication  (6): 

|5(/)|2  =  |P(/)|2-|V(/)|2  (2.2) 
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The  logarithm  now  separates  the  signal  into  two  added  terms  (6): 

log  |5(/)|J  =  log  |  P(/)|J  +  log  |V(/)|2  (2.3) 

The  power  cepstrum  of  the  input  signal  is  now  the  sum  of  the  power  cepstrum  of  the  two  convolved 
signals.  The  two  sums  may  be  easily  separated  by  filtering  if  they  occupy  different  quefrency 
ranges  (27).  Pitch  and  vocal  tract  signals  generally  do  have  different  ranges  in  quefrency  (23).  The 
pitch  of  the  voice  fluctuates  very  quickly  in  the  spectral  domain  of  speech  -  therefore  it  has  a  large 
spike  at  a  high  quefrency.  The  vocal  tract  is  made  up  of  the  lips,  mouth,  and  tongue.  These  can 
not  change  nearly  as  quickly  as  the  vocal  cords  can,  therefore  the  vocal  tract  fluctuates  slowly  in 
the  frequency  domain  and  has  a  spike  in  the  low  quefrency  region.  These  spikes  are  called  formant 
peaks  and  have  a  number  of  uses  in  speech  processing  (6).  Speech  identification  systems  do  not 
use  pitch  information;  male  and  female  speech  can  be  made  to  look  the  same  by  eliminating  the 
formant  corresponding  to  the  pitch  (2).  The  cepstrum  is  used  to  find  the  formants  of  the  pitch  so 
that  this  information  can  be  removed  from  the  speech  signal.  The  biggest  problem  with  the  power 
cepstrum  is  that  many  signals  can  not  be  directly  restored  because  the  phase  information  has  been 
lost. 

The  separated  signals  could  be  reconstructed  if  the  signal’s  phase  information  were  main¬ 
tained.  The  complex  cepstrum  accomplishes  this  by  using  the  complete  spectrum  of  the  signal 
instead  of  simply  the  power  spectrum  (10).  The  complex  cepstrum  is  defined  as  the  inverse  Fourier 
transform  of  the  logarithm  of  the  spectrum  of  a  signal  (10).  This  can  be  shown  in  equation  form 
as: 

c,  =  T-llog[S{f))  =  ±  jT  log [S(f))e»*'df  (2.4) 

The  complex  cepstrum  retains  the  phase  information  of  the  original  signals  by  taking  the  complex 
logarithm  of  the  spectrum.  The  natural  logarithm  of  a  complex  value  is  equivalent  to  the  natural 
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logarithm  of  the  magnitude  plus  the  phase  of  the  signal  (22:769)  as  shown  below: 


ln[X(/)]  =ln[|X(/)|e>Z*</>] 

(2.5) 

=  In  \X  (/)|  +  jLX(f) 

(2.6) 

The  complex  cepstrum  is  a  complex  function,  as  the  name  implies.  However,  all  complex  cepstral 
coefficients  will  be  real,  if  the  original  signal  is  real  (22).  The  real  coefficients  result  due  to  the  fact 
that  a  real  signal’s  spectrum  will  have  phase  with  odd  symmetry  and  the  resulting  inverse  Fourier 
transform  will  yield  real  coefficients.  Because  it  keeps  all  phase  and  magnitude  information,  the 
complex  cepstrum  is  reversible  -  it  is  possible  to  return  to  the  original  signal  after  performing 
filtering  operations.  Signals  may  not  only  be  separated,  like  with  the  power  cepstrum,  but  the 
signal  may  also  be  fully  reconstructed.  Parts  of  the  signal  may  be  completely  removed  with  the 
complex  cepstrum.  This  separation  and  removal  is  called  deconvolution.  A  transmitted  signed  may 
be  separated  into  its  original  signal  plus  effects  of  transmission.  These  effects  may  then  be  filtered 
and  removed  and  the  original  signal  restored  (27). 

The  complex  cepstrum  is  complicated  by  ambiguities  in  the  phase  angle,  LX(e^w)  (22:775). 
At  any  given  value  of  w,  integer  multiples  of  2x  +  w  will  yield  the  same  phase  value.  This  leads 
to  the  complex  logarithm  being  multi-valued.  The  phase  of  the  signal  will  have  discontinuities 
because  each  point  can  take  on  any  value  modulo  2x  of  its  principle  value.  There  have  been  many 
algorithms  defined  to  solve  this  phase  unwrapping  problem  (6): 

1.  A  correction  sequence  method,  where  a  correction  is  calculated  and  added  to  the  modulo  2x 
phase. 

2.  Integrating  the  phase  derivative  . 

3.  An  adaptive  numerical  integration  procedure. 

4.  A  recursive  procedure  to  remove  the  linear  phase. 
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Note  that  for  minimum  phase  signals  (no  poles  or  zeros  outside  the  unit  circle),  phase  unwrapping 
is  not  necessary.  In  this  case  the  complex  cepstrum  will  be  the  same  as  the  power  cepstrum. 

The  problems  encountered  and  the  additional  computational  complexity  are  a  significant 
drawback  to  the  use  of  the  complex  cepstrum.  Because  of  the  laborious  mathematics  required  for 
these  phase  unwrapping  techniques,  there  have  been  a  number  of  alternative  forms  for  calculating 
the  complex  cepstrum  without  dealing  with  the  phase  information.  The  differential  cepstrum  as 
described  by  Polydoros  and  Fam  (24)  takes  the  Fourier  transform,  then  finds  a  normalized  deriva¬ 
tive  of  the  transform  coefficients,  takes  the  inverse  Fourier,  and  finally  the  logarithm  to  get  the 
differential  cepstral  coefficients.  In  another  approach,  Bednar  and  Watt  estimate  the  derivative  of 
the  log-amplitude  as  well  as  the  phase  spectrum  and  then  perform  a  frequency  domain  integration 
by  dividing  in  time  (3).  This  is  similar  to  the  integration  method  of  Oppenheim  and  Schaefer  (22), 
however  Bednar  and  Watt  avoid  any  numerical  integration  or  phase  unwrapping  to  obtain  i(n) 
from  x(n).  Cepstral  coefficients  have  also  been  calculated  entirely  in  the  time  domain  in  order  to 
avoid  the  phase,  windowing,  and  aliasing  problems  of  the  Fourier  transform  (30).  This  time  domain 
cepstral  transform  is  actually  more  computationally  intensive  than  the  Fourier  based  cepstrum,  but 
it  can  yield  better  results  (30).  It  has  also  been  shown  that  homomorphic  filtering  can  be  accom¬ 
plished  without  directly  calculating  cepstral  coefficients  (11).  Khare  and  Yoshikawa  demonstrated 
a  relationship  between  the  moments  of  a  signed  sequence  and  of  its  corresponding  cepstral  sequence. 

Transforms 

Transforms  are  composed  of  matrices,  A,  which  calculate  coefficients  from  an  input  signal. 
These  coefficients  can  be  easily  processed  and  transmitted.  At  a  receiver,  the  signed  is  reconstructed 
by  performing  the  inverse  transform,  A-1,  and  putting  the  blocks  back  together  (16).  One  of  the 
most  common  transforms  is  the  block  treuisform.  The  block  transform  simply  divides  a  signed  into 
equal  blocks  and  transforms  each  block  into  a  group  of  coefficients.  Eetch  block  of  length  N  is 
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changed  by  the  transform  matrix,  A  (NxN  matrix)  so  that: 


x(n)  =  [z(nJV),  x(nN  —  1) ...  x(nN  -  N  +  1)]T  (2.7) 

This  transformation  can  also  be  written  as  a  general  linear  equation:  X  =  ATx  (16).  The  matrix 
X  is  called  the  transform  of  the  given  signal  x.  One  goal,  when  using  transform  coding,  is  to  reduce 
the  number  of  bits  per  second  required  to  transmit  a  given  signal.  This  is  accomplished  through 
careful  selection  of  A.  When  the  transform  matrix  is  orthogonal,  the  transpose  of  A  equals  the 
inverse  of  A.  We  can  then  easily  find  the  given  signal  by  the  relation:  x  =  AX  (16).  The  columns 
of  the  transform  matrix,  A,  are  called  the  basis  functions  (or  basis  vectors)  because  they  are  used 
to  transform  one  element  of  x  to  a  corresponding  element  of  X.  These  basis  functions,  a„,  are  what 
differentiate  block  transforms  from  each  other.  Orthogonal  transform  matrices  provide  us  with 
the  ability  to  easily  find  the  forward  and  inverse  transform,  because  we  do  not  need  to  actually 
calculate  the  inverse  of  the  matrix  (16). 

The  discrete  Fourier  transform  (DFT)  and  the  discrete  cosine  transform  (DCT)  are  two  of 
the  best  known  and  most  widely  used  block  transforms.  The  DFT  has  basis  functions  that  are 
complex  sinusoids  (28): 


where  N  is  the  block  size.  Like  the  Fourier  series,  the  DFT  represents  the  signal,  x[n],  as  a 
combination  of  harmonically-related  sinusoids  (16).  The  discrete  cosine  transform  (DCT)  provides 
N  different  frequencies  between  0  and  x,  where  the  DFT  has  only  N/2  frequencies.  The  DCT  basis 
is  defined  as  (28): 

a",=c(l)/l'c<,'[("+0  £]  <2'9> 


2-6 


where  c(k)  is  when  k  =  0  and  is  equal  to  one  otherwise.  The  DCT  is  so  common  in  signal 
processing  that  it  has  become  the  standard  for  the  video-phone  and  other  image  coding  applica¬ 
tions  (28). 

Performance 

Block  transforms  encounter  a  problem  called  the  blocking  effect,  which  occurs  when  signals 
are  reconstructed.  Blocking  effects  arise  due  to  the  lack  of  overlap  between  blocks  of  the  signal. 
These  effects  are  noticed  as  discontinuities  along  the  boundaries  of  an  image  or  as  extraneous  tones 
in  a  speech  signal.  These  effects  become  even  more  pronounced  at  low  data  rates  (37).  The  DFT 
and  DCT  lead  to  blocking  effects  because  of  the  independent  processing  of  each  block.  The  final 
segments  of  a  prior  block  will  generally  not  match  with  the  first  samples  of  the  next  block,  because 
the  basis  functions  (as  defined  above)  vary  from  one  block  to  the  next  (16).  Two  methods  have 
been  proposed  to  reduce  blocking  effects:  overlapping  and  filtering.  Overlapping  consists  of  blocks 
of  samples  which  simply  overlap  some  predetermined  amount.  In  this  way,  the  information  around 
the  boundaries  is  transmitted  twice  and  the  receiver  averages  the  data  in  the  overlapped  areas. 
Unfortunately,  overlapping  increases  the  bit  rate  because  extra  samples  must  be  sent.  For  filtering, 
a  low  pass  filter  is  added  at  the  receiver  to  filter  the  boundary  pixels.  Although  this  does  not 
increase  the  bit  rate,  it  tends  to  blur  the  signal  across  the  boundaries  (15). 

Malvar  explained  that  performance  of  a  multirate  filter  bank  (a  transform  coding  application) 
is  directly  related  to;  “perfect  reconstruction,  high  stopband  attenuation,  narrow  transition  width, 
high  coding  gain,  absence  of  blocking  effects,  and  fast  algorithms”  (15).  A  compromise  must  be 
reached  between  these  conflicting  goals.  Two  methods  of  coding  are  often  used:  transform  and 
subband  coding.  Subband  coding  yields  filters  of  high  performance  (ie  high  stopband  attenuation 
and  narrow  transition  widths)  but  they  are  also  of  high  complexity  (longer,  slower  algorithms). 
Transform  coding  yields  a  faster  process  but  poorer  filters,  which  lead  to  the  blocking  effects. 
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The  Lapped  Transform  (Malvar  Wavelet) 

The  lapped  transform  (or  Malvar  wavelet)  provides  an  excellent  tradeoff  between  complexity 
and  filter  performance  (15).  The  lapped  transform  has  the  same  benefits  as  the  overlapping  method, 
but  without  the  increase  in  the  bit  rate  (18).  The  LT  eliminates  the  blocking  effects  by  overlapping 
adjacent  blocks  of  the  signal,  hence  the  term  “lapped”  transform. 

First  let  the  incoming  signal,  x[n],  be  composed  of  MN  samples,  where  M  is  the  block  size. 
Traditional  overlapping  would  transform  this  signal  into  N  blocks  each  of  length  M.  The  lapped 
transform  will  instead  transform  the  signal  into  N  blocks  of  length  L,  where  L  >  M,  so  that  adjoining 
blocks  will  overlap  by  L  -  M  samples  (note  that  traditionally  L  <  2M).  The  basis  functions  have 
been  made  longer  than  the  block  length,  which  allows  for  a  smoother  transition  at  the  ends  of  the 
blocks  (see  Figure  2.2)  (16). 


ai~€j  ai  ai  +  €i  aj+i  “  ej'+i  a>+ 1  aj+i+ei+l 

Figure  2.2.  Overlapping  of  Basis  Functions 

This  lapped  transform  method  closely  resembles  the  basic  method  of  overlapping.  The  fun¬ 
damental  difference  between  the  two  is  that  the  LT  maps  the  L  samples  of  each  block  back  into  M 
transform  coefficients.  Therefore,  there  is  no  increase  in  the  data  rate  when  using  the  LT,  because 
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the  number  of  coefficients  are  now  equal  to  the  original  block  size  again  (18).  The  LT  can  also  be 
designed  as  a  perfect  reconstruction  process  -  -  if  all  of  the  transform  coefficients  are  used  at  the 
receiver,  the  reconstructed  signal  will  be  an  exact  copy  of  the  original  signal  (with  the  possibility 
of  a  pure  delay)  (37).  In  order  to  achieve  perfect  reconstruction,  the  transform  matrices  must  obey 
the  following  properties  (37): 

AAt  =  /  and  ATA  =  1  (2.10) 

where  I  is  the  identity  matrix.  This  is  the  same  as  saying  that  the  basis  functions  in  A  must  be 
orthogonal  to  each  other.  Because  the  basis  functions  overlap,  they  must  also  be  orthogonal  to  the 
basis  functions  in  the  neighboring  blocks.  The  LT  is  often  called  the  lapped  orthogonal  transform 
(LOT),  due  to  this  orthogonality  requirement  (16). 

The  LOT  is  a  better  alternative  than  the  conventional  block  transforms  for  signal  coding 
applications  because  it  eliminates  the  blocking  effects  without  increasing  the  data  rate.  The  price 
paid  is  the  increase  in  computational  complexity  of  the  transform  algorithm  (1).  Transform  coding 
generally  employs  the  discrete  cosine  transform  (DCT)  because  of  its  close  approximation  to  the 
statistically  optimal  Karhunen-Loeve  transform  (14).  The  LOT  actually  may  require  as  much  as  30 
percent  more  computations  (mainly  additions)  than  the  DCT.  However,  the  LOT  leads  to  slightly 
smaller  reconstruction  errors  than  does  the  DCT  (18).  There  have  been  a  number  of  variations  and 
extensions  to  the  theory  of  lapped  transforms  since  the  introduction  of  the  first  LOT  by  Cassereau. 
Some  of  the  other  types  of  LOTs  which  have  been  characterized  include:  (i)  Modulated  Lapped 
Transform  (MLT),  which  is  based  on  Quadrature  Mirror  Filter  Banks  (14);  (ii)  Extended  Lapped 
Transform  (ELT),  which  is  an  MLT  utilizing  a  larger  basis  function  (15);  (iii)  The  LOT  then  can 
be  generalized  with  any  orthogonal  bases  desired  (1).  The  LOT,  although  not  a  wavelet  in  the 
traditional  sense,  compliments  wavelet  packet  analysis,  as  such,  it  has  been  termed  the  Malvar 
wavelet  by  Yves  Meyer  (19). 
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The  Malvar  Wavelet  as  a  Mvttirate  FtHerbank 


Transforms  can  also  be  represented  as  a  series  of  filterbanks  (37:113).  “A  digital  filter  bank  is 
a  collection  of  digital  filters  with  a  common  input”  (37:113).  The  Fourier  transform  can  be  shown 
as  in  figure  2.3.  The  analysis  bank  splits  the  signal,  z[n],  into  M  subband  signals,  x*[n],  where 
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Figure  2.3.  The  DFT  Represented  in  Multirate  Form 

0  <  k  <  M  —  1.  The  synthesis  filters  then  recombine  the  M  subband  signals  into  a  single  signal, 
*[n].  The  Discrete  Fourier  Transform  (DFT)  is  defined  as  (37:794): 

*(*)  =  E  *WWj?  (2-11) 

n=0 

where  Wm  —  .  The  analysis  and  synthesis  banks  of  the  DFT  are  therefore  constructed  of 

the  matrix  W  which  has  MxM  elements. 

The  Malvar  wavelet  may  be  represented  as  a  multirate  filterbank  where  the  matrix  T  is  now 
MxL  (rather  than  MxM,  where  L  >  M)  as  shown  in  figure  2.4  (37:325).  This  means  that  there 
are  more  input  subbands  to  the  filter  T  than  output  subbands.  The  Malvar  wavelet,  as  stated 
previously,  will  take  the  L  samples  of  each  block  and  map  them  onto  M  transform  coefficients. 
The  |M  (in  figure  2.4)  represents  a  decimator.  A  decimator  (or  downsampler)  retains  only  those 
samples  of  z[n]  which  occur  at  time  equal  to  multiples  of  M  (37:100).  For  example,  if  M=2  then 
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Figure  2.4.  The  Malvar  Wavelet  Represented  in  Polyphase  Form 


the  signal  would  be  broken  into  its  even  and  odd  components.  The  |  M  then  puts  the  samples  back 
into  their  proper  place  in  the  data  stream  when  the  signal  is  reconstructed  at  the  receiver. 

By  downsampling  prior  to  transmission,  the  computation  rate  through  a  communications 
satellite  can  be  reduced.  Therefore  slower,  cheaper,  more  efficient,  and  more  reliable  elements  can 
be  used  on-board  the  satellite  (32).  This,  in  turn,  would  mean  that  the  lose  of  an  element  would 
have  less  of  an  effect  on  the  mission. 

Generalized  Malvar  Wavelet 

Windowing,  often  employed  in  filtering  and  spectral  estimation,  can  result  in  irreparable 
distortions  in  the  transformed  signal.  By  utilizing  the  Malvar  wavelet,  any  signal  distortion  resulting 
from  the  transformation  can  be  eliminated  at  reconstruction.  Suter  and  Oxley  have  presented  a 
Generalized  Malvar  wavelet  (GMW)  which  places  conditions  on  the  windows  and  orthonormal  basis 
and  allows  for  perfect  reconstruction  (34).  The  basis  is  defined  by  performing  an  odd  and  even 
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extension  of  the  orthogonal  basis  function,  />,*(*),  about  the  end  points  of  the  window  such  that: 


/;,*(*)  =  { 


0  ,  —oo  <  x  <aj  —  tj 

~fjA2ai  ~  x)  .  ai  ~  *i  <  x  <  a, 

fj,k(,x )  i  aj  ^  x  5:  “j+1 

fi,k(2ai+ 1  —  z)  •  ai+i  <  *  <  ai+i  +  ei+i 

0  ,  aJ+i  +  fj+i  <  x  <  oo. 


(2.12) 


where  aj  and  aj+i  are  the  end  points  of  the  frame  and  tj  and  e;-+i  are  the  overlaps  about  the 
two  points.  An  associated  weight  function,  Pj(x),  is  then  used  to  ensure  perfect  reconstruction  is 
possible.  All  of  the  basis  functions  investigated  in  this  thesis  do  not  require  a  weight  function;  let 
Pj(x )  =  1.  The  normalized  window  of  the  GMW  is  chosen  such  that  (8): 


Wj(x)  =  1 

for 

a,  -f  Cj  <  x  <  Oj+i  -  c;+1 

o 

II 

s' 

for 

0;  -  €;•  >  x  >  Oj+l  + 

Wj(aj  ~  a)-  +  <r) 

for 

<r  £[-(}, (j] 

w](x)  +  U)J_x(x)  =  1 

for 

aj  —  Cj  <  x  <  aj  +  tj 

(2.13) 

The  Generalized  Malvar  wavelet  (generalized  lapped  orthogonal  transform)  is  thus  a  lapped 
transform  which  allows  for  many  variations  within  the  same  transform.  Most  transforms,  like  the 
DFT  or  DCT,  maintain  the  same  block  length  and  basis  function  throughout  the  data  set.  Suter 
and  Oxley  derived  the  GMW  based  on  being  able  to  vary  the: 

1.  Lengths  of  the  Windows. 

2.  Orthonormal  Basis  Functions  (A). 

3.  Overlap  Between  Adjacent  Windows. 

4.  Window  Amplitudes. 
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The  goal  of  the  generalized  Malvar  wavelet,  like  all  transform  coding,  is  to  determine  a  set  of 
coefficients  which  may  then  be  used  to  reconstruct  the  original  data.  The  benefit  of  being  able  to 
vary  the  parts  of  the  transform  is  that  the  transform  may  be  tailored  to  the  data  set.  For  example, 
the  number  of  coefficients  required  to  model  a  slowly  varying  signal  can  be  reduced  by  increasing 
the  length  of  the  window  (32). 

Two  Malvar  wavelets  were  studied  for  this  thesis  effort.  A  Complex  Valued  Malvar  wavelet, 
which  is  very  similar  to  the  short  time  Fourier  transform  because  it  utilizes  an  exponential  basis 
function  (32): 

f(x)  =  e>2*mx  (2.14) 

With  the  weight  function  ,  p(x),  equal  to  one.  The  Complex  Malvar  wavelet  also  keeps  a  constant 
length  window  throughout  the  data  set  as  well  as  constant  amounts  of  overlap  between  the  windows. 
The  second  transform  is  a  Generalized  Real  Valued  Malvar  wavelet.  This  wavelet  uses  a  sinusoidal 
basis  function: 

/*(*)  =  (*(*  +  V2)  [£])  (215) 

with  p(x)  =  1.  This  transform  is  a  generalized  wavelet  because  it  may  use  variable  length  windows, 
and  variable  amounts  of  overlap  within  any  given  data  set  (34). 
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III.  Complex-Valued  Malvar  Wavelets 


Introduction 

The  Balian-Low  Theorem  states  that,  if  a  complex-valued  exponential  basis  is  windowed, 
the  resulting  transform  will  not  be  localized  in  time  or  not  be  localized  in  frequency.  As  a  result, 
orthogonality  of  the  transform  can  not  be  maintained.  This  represents  a  significant  constraint 
to  signal  processing  because  of  the  central  role  of  the  Short  Time  Fourier  Transform.  Although 
this  represents  a  “no-go”  condition,  by  changing  the  underlying  assumptions  a  positive  result  is 
possible.  Namely,  if  the  functions  have  extensions  outside  the  interval  and  windows  are  defined 
with  constraints,  an  orthogonal  windowed  transform  is  possible. 

Suter  and  Oxley  have  formulated  an  orthogonal  Malvar  wavelet  that  has  a  complex-valued 
exponential  basis  (33).  Moreover,  the  resulting  algorithm  has  the  same  “Big  Oh”  complexity  as 
the  conventional  Fast  Fourier  Transform  (FFT). 

Complex-  Valued  Malvar  Wavelet  Algorithm 

Suter  and  Oxley  (33)  generate  the  complex-valued  Malvar  wavelet  (CMW)  with  eight  primary 
steps  which  are  summarized  below: 

(1)  Define  a  Basis:  An  exponential-like  basis  can  be  applied  to  the  Malvar  wavelet  and 
complete  orthogonality  will  be  preserved,  even  after  windowing.  The  basis  function  is  first  defined 
as  the  even  and  odd  extensions  of  an  exponential  function: 

0  ,  — oo  <  x  <  n  —  e 

e-}2rmx  J  n  —  €  <  X  <Tl 

fn,m(x)  =  e*2™*  ,  n  <  x  <  (n  + 1)  (3-l) 

t  (B  +  l)<*<(n+l)  +  e 

0  ,  (n  +  1)  +  e  <  *  <  oo. 
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Where  [n,  n  +  1]  is  the  nth  interval,  and  e  is  the  amount  of  overlap  between  windows.  The  basis 
function  definition  is  equivalent  to  even  and  odd  extending  both  the  real  and  imaginary  parts  of  the 
function  beyond  (n)  and  (n+1).  The  basis  functions  and  their  extensions  are  shown  in  figures  3.1a 
and  b.  Note  that  the  values  of  n  and  n+1,  in  the  figure  and  the  calculations,  have  been  normalized 
to  zero  and  one. 

(2)  Define  a  Window:  The  transform  will  have  a  single  window  function  that  will  be  used 
over  each  interval.  The  window,  g(x)  (with  n  =  0),  is  defined  as: 


9(x) 


0  ,  — oo  <  z  <  0  —  t 

sin(^{x  +  e})  ,  0-c<x<0+f 

'  1  ,  0  +  f<X<(0+l)-£ 

cos(£{x  -  (1)  +  £})  ,  (0  +  1)  —  f  <  X  <  (0+  l,  +  ( 

0  ,  (0  +  1)  +  f  <  x  <  oo. 


(3.2) 


Now  the  basis  is  equal  to  un  m  =  fn,m(x)g(x  -  n).  This  is  an  orthonormal  basis  set  for  L2(S). 

(3)  Multiply  By  the  Window  and  Define  the  Coefficients,  a:  The  complex  Malvar  wavelet 
coefficients  are  calculated  by: 


a 


n,m 


-L 

-L 


{ st  un,m) 
n+l+e 


s(x)g(x  -  n)/n>m(x)dx 


n-< 

n+1 


Mx)eJ2Tm*dx 


(3.3) 
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where  s(x)  is  the  input  signal  and  g(x)  and  f{x)  are  as  defined  above  and  hn  is  defined  as: 


Is(x)g(x  —  n)  +  s(2n  -  x)g(n  -  x)  ,  n  <  x  <  n  +  e 

a(x)  ,  n  +  r<x<(n  +  l)-e  (3-4) 

a(x)g(x  -  n)  -  a(2(n  +  1)  -  x)g(n  +  2-  x)  ,  (n  +  1)  -  c  <  x  <  (n  +  1) 

(4)  Define  a  New  Array:  Let  Hnj  =  h„(n  +  61).  Where,  Hn, o  =  h„(n)  denotes  the  Folded 
Data  at  the  start  of  the  interval,  and  Hn,N  =  hn(r»  +  1)  denotes  the  Folded  Data  at  the  end  of  the 
interval.  The  trapezoidal  rule  can  then  be  used  to  approximate  the  integration: 


N-l 


«n.m  *  <m  =  \Hnfi  +  £  \ Hn,N 


(3.5) 


1=1 


From  the  hn  equations,  it  is  clear  that  Hn  pj  will  equal  zero.  Now  let: 


Pn,l  =  < 


ffn.O  ,  /  =  0 

2Hnj  ,  1</<AT-1 


(3.6) 


(5)  Perform  an  Inverse  FFT:  The  summation  from  step  4  thus  becomes  an  N  point  inverse 
Fourier  transform  of  fl.  The  spectral  coefficients,  a,  are  calculated  with  a  Fourier  transform  of  the 
windowed  and  scaled  input  signal: 


an,m 


1=0 


(3.7) 


The  Reconstruction  of  the  signal  is  accomplished  by  reversing  the  above  steps: 


(1)  Perform  a  Forward  FFT:  The  signal  may  now  be  reconstructed,  given  the  coefficients 


fin ,/  —  2  • 


(3.8) 


(2)  Scale  the  Values:  Scaled  the  0  values  back  to  Hn,L- 


Hn,i  = 


0n,O  ,  /  =  0 

(i)A.,!  .  1<1<N-1 


(3.9) 


(3)  Divide  By  the  Window:  Finally,  the  reconstructed  signal,  Sn>/,  is  calculated  by  dividing 
back  out  the  original  window,  g(x ): 


Sn,,= 


< 


g{x  -  n)Hnj  -  g(n  -  x)Hn~i,N~ I 

HnJ 

g(x  -  n)Hn>l  +  g((n  +  2)  -  x)H„+iiN-i 


,  n  <  x  <  n  +  t 
,  n  +  c<r<(n+l)  —  e 
,  (n+1)  — e<i<n+l 


(3.10) 


Software  Implementation 

It  is  important  that  follow-on  students  be  able  to  understand,  utilize,  and  modify  the  Mal- 
var  wavelet  programs.  The  Ada  programming  language  is  well  known  for  being  easily  modified 
and  maintained.  Therefore,  the  complex  Malvar  wavelet  was  implemented  using  Ada.  The  com¬ 
plex  transform  program  will  correspond  to  the  sinusoidal  Malvar  wavelet  Ada  program  written  by 
Raduenz  (26).  Thus  all  of  the  Malvar  code  will  be  in  the  same  language  and  be  easily  compared 
and  contrasted.  This  section  will  step  through  each  package  of  the  wavelet  program  to  outline  what 
it  is  doing  and  how  the  package  relates  to  the  rest  of  the  code.  A  complete  listing  of  the  source 
code  is  listed  in  Appendix  B. 


GetJ'ilc.Data:  The  Gct-File-Data  procedure  simply  reads  in  the  required  input  files  and 
then  prepares  the  data  vectors  to  be  used  by  the  rest  of  the  program.  The  user  must  give  the  name 
of  an  input  file  to  be  used  by  this  procedure.  This  input  file  must  contain  the  following: 

1.  The  filename  containing  the  data  set  to  be  processed. 

2.  The  filename  for  the  spectral  coefficient  data. 

3.  The  filename  for  the  reconstructed  data  set. 

4.  The  number  of  data  points  within  the  data  set  (or  number  of  points  wished  to  be  processed). 

5.  The  window  size  to  be  used  (this  must  be  an  integer  multiple  of  the  size  of  the  data  set). 

6.  The  amount  of  overlap  to  be  used  between  intervals. 

An  example  input  file  is  shown  in  figure  3.2.  It  reports  a  data  file  called  “input.dat”  and 
requests  the  spectral  coefficients  to  be  placed  in  “alpha.dat”  and  the  reconstructed  data  set  into 
the  file  called  “output.dat” .  The  file  is  1000  points  long  (or  just  want  to  process  1000  points  of 
a  larger  data  set).  The  window  size  is  200  points.  Therefore  there  will  be  5  windows  of  data: 
{(0,200),  (200,400),  •  •  •,  (900,1000)}.  Now  the  first  point  (0)  and  the  last  point  (1000)  are  actually 
the  same  partition  point.  There  is  no  data  which  resides  in  point  zero;  it  is  used  only  as  a  place 
holder  for  the  overlap  region.  The  first  window  will  overlap  with  the  last  window  and  the  last 
will  overlap  with  the  first.  This  overlap  has  been  accomplished  without  copying  any  data.  This 
initial  overlap  could  be  modified  so  that  the  first  window  would  overlap  with  itself,  thus  allowing 
for  a  continuous  stream  of  data.  The  program  calls  in  the  correct  data  positions  to  accomplish  the 
folding  and  unfolding  of  all  of  the  windows,  including  the  first  and  the  last.  Finally,  the  last  value 
tells  us  that  there  are  10  points  of  overlap  between  intervals.  The  partitions  are  established  after 
all  of  the  data  is  read  into  the  program.  The  program  determines  where  all  of  the  intervals  are 
going  to  start  and  stop.  These  correspond  to  the  positions;  {N(n),  N(n+1),  N(n+2),  •  •  •  }.  The 
last  thing  Gtt-FiltJ)aia  does  is  to  read  in  the  actual  data  set  into  a  vector  called  “Data(j)”. 

Multiply  Jy.  Windows:  This  procedure  performs  the  “folding”  of  the  data  by  multiplying 
the  data  by  the  window  as  defined  by  equation  3.4.  The  window  and  data  values  are  passed  to  this 
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input  .dai 
ilphiJtt 
output  .djU 
1000 
200 
10 

Figure  3.2.  Example  Input  File 

procedure  from  the  master  package,  Do.  Work.  The  procedure  steps  through  equation  3.4  starting 
with  the  left  edge  of  the  interval,  then  the  center,  and  finally  the  right  edge.  There  are  unique 
right  edge  equations  for  the  first  window.  This  is  due  to  the  fact  that  data  must  be  pulled  from 
the  other  end  of  the  data  file  (and  not  just  the  adjoining  interval).  For  the  same  reasons,  there  are 
special  equations  for  the  left  edge  when  using  the  last  interval.  The  folded  data  is  then  sent  back 
to  the  master  package,  Do.Work. 

Compute-Coefficients:  This  procedure  receives  the  folded  data  in  segments  (from  N(n)  to 
N(n+1)-1,  or  0  to  N-l).  It  first  scales  the  coefficients  as  in  equation  3.6  and  converts  the  data  into 
complex  notation.  It  then  sends  the  complex  data  to  the  FFT  package  to  perform  an  N  point  (size 
of  the  interval)  inverse  fast  Fourier  transform.  The  output  of  the  FFT  is  then  scaled  by  |  per  the 
algorithm.  This  complex  data  (form  0  to  N-l)  is  then  sent  back  to  Do.Work. 

ReconsiructionJTT:  The  reconstruction  FFT  package  simply  takes  in  the  N  points  passed 
from  the  Compute. Coefficients  and  sends  it  to  the  FFT  package.  Within  the  FFT  package,  an  N 
point  forward  FFT  is  performed.  The  very  first  coefficient  is  then  doubled  because  of  the  way  0 
was  defined  (see  equation  3.6)  and  the  entire  segment  is  sent  back  to  Do.Work. 

Divide  Jiy.Windows:  Divide  Jy. Windows  calculates  the  reconstructed  signal  as  defined  in 
equations  3. 10.  As  shown  in  the  algorithm,  this  package  calculates  the  segment  from  (n  —  c)  to 
((n  4-  1)  4-  e).  If  the  regular  interval,  from  n  to  (n+1),  were  calculated  then  coefficients  from  the 
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previous  interval  would  be  altered  prior  to  being  used  to  calculate  the  values  for  the  next  interval. 
As  in  Multiply  Jy.  Windows,  there  are  special  equations  for  the  first  window.  The  first  window  will 
be  pulling  data  from  the  end  of  the  data  file.  The  last  window  does  not  need  any  special  equations 
because  we  only  calculate  up  to  the  end  of  file  minus  the  overlap.  These  “missing”  points  (from  the 
end  of  file  -  overlap  to  the  end  of  file)  have  been  calculated  with  the  first  window  and  are  located 
from  point  (0)  to  (0  —  e). 

Error. Out:  This  procedure  calculates  the  reconstruction  errors  as  defined  later  in  this  chap¬ 
ter.  It  receives  the  reconstructed  and  original  data  sets  and  calculates  the  three  numerical  error 
criteria: 

1.  Relative  L 2  Error  based  on  equation  3.11. 

2.  Loo  Error  based  on  equation  3.12. 

3.  Root  Mean  Square  Error  based  on  equation  3.13. 

4.  Normalized  Renyi  Information  Value  based  on  equation  3.15. 

Do.  Work:  This  package  controls  the  entire  wavelet  procedure.  This  procedure  first  calcu¬ 
lates  the  window  based  upon  the  set  of  equations  3.2.  This  window  is  normalized  from  zero  to  one 
and  is  used  for  all  folding  and  unfolding  of  the  data  set.  Each  data  segment  is  processed  in  the 
main  loop  of  this  procedure.  The  main  loop  is  defined  from  (0)  to  (Partitions’last-1)  this  means 
that  the  program  will  start  with  the  first  window  (point  0  to  the  first  partition)  and  end  with  the 
last  window  (the  last  partition  -  1  to  the  last  partition,  which  is  the  end  of  the  file).  In  this  way,  the 
finite  length  signal  can  be  calculated  without  copying  data  points  or  processing  segments  multiple 
times. 

Do.  Work  calls  each  procedure  as  it  is  needed  by  the  algorithm.  This  main  procedure  sends 
all  of  the  necessary  data  points  and  receives  the  transformed  data  back  from  the  procedures.  The 
main  loop  is  started  over  prior  to  the  Divide  JBy.  Windows  procedure,  because  the  reconstruction 
of  the  first  window  of  the  data  set  requires  the  transformed  data  from  the  end  of  the  file.  After 
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all  steps  of  the  wavelet  are  complete,  Do.  Work  shifts  the  fust  c  values  of  the  reconstructed  data 
to  the  end  of  the  reconstructed  data  file.  This  is  done  because  Divide  Jdg.Windows  calculates  the 
data  based  on  shifted  windows  (n-€<x<(n  +  l)-f  instead  of  n  <  x  <  (n  +  1))-  Do.  Work  then 
prints  all  of  the  information  out  to  the  requested  data  files. 

Called  Packages  The  complex  Malvar  wavelet  program  references  multiple  outside  called 
packages,  including: 

1.  Text  Jo 

2.  ComplexJ’kg 

3.  MathJiib 

4.  VectorJ’ackage 

5.  TypeJ>ackage 

6.  Print J>ackage 

7.  FFTJ’ackage 

These  packages  are  discussed  and  documented  by  Raduenz  (26:4-10)  and  are  quite  integral 
to  the  workings  of  all  of  the  wavelet  programs. 

Error  Criteria 

When  testing  and  evaluating  speech  coding  systems,  one  is  usually  interested  in  two  measures 
-  quality  and  intelligibility.  These  two  measures  are  closely  related  but  not  identical.  If  the  quality 
of  a  signal  is  high  then  the  intelligibility  will  also  be  high.  However,  the  opposite  is  not  true,  it  is 
possible  to  have  very  intelligible  speech  that  is  not  of  high  quality.  Speech  synthesis  systems  produce 
speech  that  is  readily  understandable,  but  they  have  a  definite  synthetic  sound  which  degrades  their 
quality.  Both  quality  and  intelligibility  are  difficult  to  quantify.  Coders  are  not  normally  tested 
using  numerical  methods  but  with  a  group  of  listeners.  The  listeners  indicate  which  rhyming  word 
from  a  given  set  was  perceived.  The  number  of  words  correctly  perceived  provides  a  score  that 
indicates  the  intelligibility  of  the  system. 
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Numerical  error  criteria  were  utilized  to  determine  the  basic  capabilities  of  each  system.  These 
error  criteria  were  very  useful  when  simple  sine  waves  and  other  known  signals  were  input  into  the 
systems.  Six  different  criteria  were  used;  four  of  these  are  numerical  and  two  are  more  subjective 
in  nature. 

The  relative  L2  measure  is  the  error  of  the  L2  norm  of  a  signal  (approximated  by  n  coefficients) 
relative  to  the  L2  norm  of  the  original  signal  (using  all  coefficients)  (32).  The  Loo  criteria  is  a 
measure  of  the  point-wise  peak  error;  the  point-wise  maximum  difference  between  two  signals  (32). 
Coifman  used  the  relative  L2  error  criteria  when  he  researched  the  link  between  compression  of  a 
signal  and  the  basis-set  of  a  given  window  (7).  The  L2  and  L^  are  defined  as: 


L2  = 


actual2(x)  —  recon2 (x) \ 

actuaP(x) 


(3.11) 


Loo  =  Max  [|arfua/(z)  -  recon(z)|] 


(3.12) 


The  normalized  root  mean  square  (RMS)  criteria  is  the  basic  error  criteria  that  compares 
the  point-wise  difference  between  two  signals.  It  relates  the  total  power  difference  between  the  two 
signals.  The  normalized  RMS  is  defined  as: 


RMS  = 


l^£>dua/( 


z)  -  recon(z))2 


(3.13) 


The  final  numerical  error  criteria  used  is  the  Renyi  Information  Value  (RIV).  Information 
measures  are  normally  used  to  describe  the  distributions  of  random  variables  (40:146).  However 
this  measure  has  the  potential  to  give  some  real  insight  as  to  how  similar  two  signals  are.  The  Renyi 
value  relates  the  total  amount  of  “information”  residing  in  the  signal.  When  a  signal  is  transformed, 
compressed,  and  reconstructed,  it  should  loose  some  of  the  original  signals  information.  The  general 
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definition  of  the  ath  order  Renyi  information  value  for  a  discrete  function,  p,  is  (40:147): 


Ra{p)  =  y~  lo*2  Y  p“(x) 

x 


(3.14) 


where  0<a<lorl<a<oo.  In  order  to  obtain  a  relation  between  the  original  signals  informa¬ 
tion  and  the  reconstructed  signals  information,  the  second  order  Renyi  value  of  the  reconstructed 
signal  was  divided  by  the  second  order  Renyi  value  of  the  original  signal  (32).  Using  this  method, 
a  value  of  1.0  would  represent  perfect  reconstruction  (with  respect  to  the  amount  of  information 
present)  and  0.0  would  represent  a  total  loss  of  the  signal  (100%  compression).  The  Renyi  measure 
will  take  the  form: 


R2(actual,  recon) 


[logs  E«  recon2(a:) 
log2  Ey  actually) 


(3.15) 


The  final  two  error  tests  are  a  more  subjective  criteria,  but  may  have  been  the  most  important 
tests  of  all.  First,  the  reconstructed  signals  were  plotted  and  visually  compared  to  the  original 
signal.  The  plots  gave  a  great  amount  of  insight  into  the  amount  of  noise  in  the  reconstructed 
signal  and  how  much  information  was  lost.  Then,  in  the  case  of  speech  signals,  the  reconstructed 
signal  was  played  back  via  Entropies  Waves-f  •  This  is  the  ultimate  test  for  a  digitized  speech  signal. 
A  signal  earn  have  the  lowest  reconstruction  error  possible,  but  if  it  does  not  sound  correct,  then 
the  reconstruction  has  been  poor. 


Validation 

The  program  was  debugged  with  two  simple  input  files;  a  unit  step  and  a  sine  wave.  Once 
these  were  reconstructed  correctly,  the  size  of  the  data  file,  size  of  the  windows,  and  the  amount  of 
overlap  were  varied  to  ensure  that  all  of  these  combinations  yielded  perfect  reconstruction.  More 
complex  signals,  like  speech,  were  used  as  the  final  test  of  the  program’s  reconstruction  abilities. 
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Table  3.1  presents  the  reconstruction  errors  associated  with  the  three  signals;  unit  step, 
sinusoid,  and  speech.  These  signals  were  broken-down  and  reconstructed  using  a  sample  size  of 
2000,  with  250  point  windows  and  10  points  of  overlap  between  windows.  The  last  entry  on  the 
table  is  for  an  entire  sentence  of  speech  (28250  points  at  8KHz),  with  250  point  windows  and  10 
point  overlap.  The  table  shows  that  each  of  the  four  signals  were  reconstructed  to  within  the 
machine  round-off  of  the  original  data  set.  Note  that  the  RIV  error  is  plotted  as  1.0  plus  the  given 
value  in  the  table.  This  was  done  because  all  of  the  RIV  outputs  were  0.999-  •  -999,  which  does  not 
help  to  distinguish  the  errors.  The  error  rises  slight  as  the  signals  get  more  complicated  and  as 
the  amplitudes  increase.  The  worst  RMS  error  was  the  utterance  of  the  word  “dark.”  This  is  a 
relatively  short,  high  amplitude,  and  quickly  varying  signal  that  was  still  reconstructed  perfectly 
(RMS  error  =  4.216a:10-10).  The  errors  introduced  when  reconstructing  an  entire  sentence  were 
lower  because  the  numerical  error  was  spread  out  and  averaged  across  a  much  larger  number  of 
samples. 


Input  Data 

Relative  Li 

£oo 

RMS 

WJSEtfEXSM 

Unit  Step 

Sinusoid 

Utterance  of  “DARK” 
Speech  Sentence 

in 

mi 

H3 

gjf 

Table  3.1.  Output  Errors  of  the  CMW  Program  with  Different  Input  Signals 


Table  3.2  shows  how  the  error  changes  with  increasing  overlap.  The  signal  used  was  the 
2000  point  unit  step  with  250  point  windows.  Early  runs  showed  that  the  majority  of  the  error 
is  found  in  the  overlap  regions.  This  overlap  error  is  provided  in  table  3.3,  which  shows  the  given 
reconstruction  around  an  end  point  (N)  of  a  unit  step  function.  Thus,  for  signals  that  are  not 
compressed,  as  the  amount  of  overlap  increases  so  will  the  calculation  error. 

The  Malvar  wavelet  program  not  only  provide  very  low  numerical  errors  but  also  the  graphic 
displays  and  the  audio  reconstructions  were  perfect.  The  reconstruction  of  the  unit  step  and  the 
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Amount  Overlap 

Relative  L i 

Lqq 

RMS 

mnwfwiwmu 

1  point 

TiraiiicFs- 

2.442*10-m 

1.807*10~l® 

0.000*10-“ 

5% 

2.482*10"® 

1.908*  10- 10 

5.628*10~13 

— 0.081*10-12 

10% 

3.728*10"® 

1.984*10-10 

8.446*  10- 13 

— 0.183*10"12 

20% 

5.367*10-® 

2.017*10-1° 

1.215*10~12 

— 0.379*10-12 

30% 

6.611*10-® 

2.028*10-10 

1.497*10~12 

— 0.575*10-12 

40% 

7.656*10"® 

2.034*10-1° 

1.734*10~12 

— 0.771*10-12 

50% 

8.575*10"® 

2.037*10"10 

1.942*10~12 

— 0.968*10-12 

Table  3.2. 


Output  Errors  of  the  CMW  Program  with  the  Unit  Step  and  Varying  Amount  of 
Overlap 


Position 

Reconstructed  Value 

N  -  6 

9.99999999999999E-01 

N  -  5 

1.00000000000002E+00 

N  -  4 

9.99999999997335E-01 

N  -  3 

9.99999999995931E-01 

N  -  2 

9.99999999995932E-01 

N  -  1 

9.99999999997335E-01 

N 

1.00000000000002E+00 

N  +  1 

1.00000000001683E+00 

N  +  2 

1 .0000000000 1252E+00 

N  +  3 

1.00000000000798E+00 

N  +  4 

1.00000000000367E+00 

N  +  5 

1.00000000000001E+00 

N  +  6 

1.00000000000000E+00 

Table  3.3.  Reconstruction  of  the  Unit  Step  Function 


utterance  of  the  word  “dark”  are  shown  in  figures  3.3  and  3.4.  There  were  no  distortions  of  the 
signal,  even  when  the  original  and  reproduced  signal  were  compared  in  very  small  segments.  The 
ASCII  data  files  may  be  transformed  into  speech  files  (.sd)  to  be  played  on  Entropies  Waves-t-  using 
the  following  commands: 

1.  addfea  -f  frequency .plot  -t  FLOAT  -c  comment  filename. asc  filename. out 

2.  tofeasd  -f  frequency .plot  -R  8000  filename. out  filename. sd 

When  played  through  the  Ariel  audio  interface,  the  original  and  recreated  words  and  sentences 
could  not  be  differentiated  from  each  other. 
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Figure  3.4.  (a)  Utterance  of  the  Word  “Dark”  and  (b)  The  Complex  Malvar  Wavelet 

Reconstruction 
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Compression  Methods 


There  are  many  different  methods  which  can  be  used  to  compress  the  number  of  coefficients 
used  for  the  transmission  of  a  signal.  Switzer  studied  eight  different  ways  to  select  the  most  critical 
values  of  a  transformed  signal  (35:31-39).  The  first  scheme  was  a  simple  max-picking  method,  where 
the  M  largest  values  were  used.  The  second  scheme  was  a  max-picking  method  with  the  frequencies 
above  600  Hz  being  boosted.  The  third  selection  scheme  used  a  peak-picking  method,  where  the 
max  values  are  selected  and  the  12  nearest  neighbors  were  then  eliminated  (to  allow  for  more  local 
maxima  to  1  e  selected).  The  fourth  used  a  max- peak  selection  with  the  four  nearest  neighbors 
eliminated.  The  other  methods  used  these  plus  a  method  for  neighbor  estimation.  He  found  that 
the  basic  max-picking  technique  performed  better  than  most  of  the  other  methods  and  worse  than 
just  two  of  the  more  complicated  methods  (35:70).  Raduenz  also  performed  similar  experiments 
with  the  LOT  (Malvar  wavelet)  and  came  to  the  same  conclusion  that  a  basic  max-picking  method 
performs  just  as  well  as  most  more  elaborate  methods  (26:4-15). 

Absolute  Thresholding:  Selecting  the  maximum  values,  as  the  max-picking  methods  of 
Switzer  and  Raduenz  do,  is  the  same  as  eliminating  all  of  the  transform  coefficients  under  a  given 
threshold.  The  Lapped  Transform  program  was  originally  written  so  that  a  threshold  could  be 
arbitrarily  selected.  Any  coefficients  under  this  threshold  would  be  set  to  zero  and  not  used  in  the 
reconstruction  of  the  signal.  For  the  experiments  being  performed,  a  given  PERCENT  compres¬ 
sion  had  to  be  obtained,  so  the  program  was  modified  to  allow  the  user  to  select  a  given  percent 
compression.  If  one  requires  85%  compression  (85  out  of  every  100  coefficients  are  set  to  zero), 
the  code  will  recursively  and  iteratively  determine  the  amount  of  thresholding  required  to  achieve 
85%  compression.  The  source  code  listed  in  Appendix  E  was  added  to  the  Do  Work  routine  of  the 
Malvar  wavelet  program. 

The  thresholding  package  first  determines  the  number  of  zeros  required  to  obtain  the  required 
percent  compression.  It  then  selects  an  initial  threshold  equal  to  the  percent  compression  times 
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the  maximum  coefficient  value  (so  that  the  initial  threshold  increases  as  the  compression  ratio 
increases).  The  procedure  then  thresholds  the  coefficients  by  setting  a  temporary  value  to  zero  and 
compares  the  number  of  zeroed  temporary  values  to  the  number  of  zeros  required.  If  there  are  too 
many  zeros,  the  threshold  is  decreased  and  if  there  are  too  few  zeros,  the  threshold  is  increased. 
Once  the  correct  number  of  zeros  are  achieved,  the  program  drops  out  of  the  infinite  loop  and 
the  coefficients  are  multiplied  by  the  temporary  values  (1.0  if  value  is  retained  or  0.0  if  value  is 
thresholded). 

Data  Rates 

The  overall  goal  for  coding  schemes  is  to  reduce  the  data  rates  necessary  to  transmit  essentially 
the  same  information.  How  is  the  compression  ratio  related  to  the  actual  data  rate?  There  are 
two  areas  that  need  to  be  addressed  in  order  to  calculate  the  data  rate.  The  first  is  how  many 
bits  will  be  required  to  represent  the  spectral  coefficients  that  will  be  transmitted.  The  second 
factor  is  determining  the  overhead  required  to  relate  the  coefficient  to  its  position  in  the  window 
(for  example,  was  the  sample  the  first,  last  or  some  other  harmonic  within  a  given  frame).  The 
typical  sampling  rate  for  communications  equipment  is  at  8  kHz  (2),  therefore,  all  experiments  were 
performed  at  the  8  kHz  rate.  The  DARPA  TIMIT  Phonetic  Speech  database  used  for  this  research 
is  sampled  at  16  kHz  so  it  was  down  sampled  to  8  kHz  prior  to  processing. 

The  coefficients  of  the  Malvar  wavelet  have  a  symmetry  such  that  only  half  of  the  coefficients 
need  to  be  used.  For  a  window  of  size  N,  the  complex-valued  Malvar  wavelet,  like  the  FFT,  will 
yield  N/2  unique  magnitude  and  N/2  unique  phase  coefficients.  The  amplitude  and  position  of 
these  coefficients  must  now  be  transmitted  using  the  least  number  of  bits  possible. 

Amplitude  Information:  Initial  experiments  were  performed  with  full  32  bit  accuracy  of 
the  coefficients.  Compression  of  90%  to  95%  of  the  coefficients  (magnitude  and  phase)  resulted  in 
speech  files  that  had  a  mechanical  tone  of  voice  with  ringing  bell-like  artifacts.  The  speech  may 
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have  warbled  and  had  a  “spacey”  sound  to  it,  but  it  was  still  intelligible.  The  higher  compression 
was  slightly  lower  quality  with  more  sudden  starts  and  stops  and  more  evident  clicking,  but  the 
overall  intelligibility  was  the  same.  The  number  of  bits  to  code  the  magnitude  and  phase  ampli¬ 
tudes  were  then  reduced  and  the  effects  on  the  numerical  and  graphical  errors  and  on  the  speech 
quality  /intelligibility  were  recorded.  Lloyd’s  algorithm  was  used  to  determine  the  quantization 
levels  of  the  magnitude  and  phase  (13). 

Lloyd’s  algorithm  starts  by  computing  the  mean  of  all  the  elements  in  the  data  to  form  a 
codebook  of  length  one.  Larger  codebooks  are  then  generated  by  splitting  the  previous  codebook 
(replacing  each  codeword  by  two  codewords,  separated  by  a  small  number),  and  going  through  a 
process  of  improving  the  codebook  until  convergence  is  reached.  This  process  works  as  follows:  The 
squared  distance  is  computed  between  every  element  of  data  and  every  codeword.  Each  codeword 
is  then  replaced  by  the  mean  of  all  the  elements  of  data  for  which  that  codeword  is  the  nearest 
neighbor.  The  mean  square  error  (MSE)  that  results  when  each  element  of  data  is  replaced  by  the 
nearest  element  in  the  codebook  is  computed.  The  improvement  loop  stops  when  the  difference 
between  the  old  distortion  (MSE)  and  the  new  distortion,  divided  by  the  old  distortion,  becomes 
less  than  the  convergence  threshold. 

Optimally  quantizing  a  set  of  coefficients  and  then  using  the  same  data  to  transform  the 
speech  file  would  not  be  a  fair  test  because  only  one  set  of  quantization  levels  is  used  in  a  voice 
coder  for  ALL  speakers.  To  avoid  the  use  of  an  optimal  quantization,  Lloyd’s  algorithm  was  used 
to  quantize  one  speech  data  file  and  those  levels  were  used  to  transform  a  second  speech  file.  In  the 
first  example,  sixteen  levels  (4  bits)  were  used  for  both  the  magnitude  and  the  phase  information. 
For  4  bit  quantization  of  magnitude  and  phase  information,  the  cutoffs  are  more  dramatic  at  the 
start  and  stops  than  the  32  bit  quantization,  however  the  overall  quality  is  almost  as  good  as  the  32 
bit  reconstruction.  Experiments  involving  several  other  combinations  of  magnitude  and  phase  bits 
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(denoted  magnitude/ phase  bits)  were  performed.  The  errors  associated  with  a  few  representative 
combinations  of  magnitude  and  phase  bits  are  summarized  in  table  3.4. 


Mag/Phast  Bits  | 

Relative  Li  @90% 

m.wmi 

Relative  Li  @95% 

RMS  @95% 

0.21783 

0.47239 

0.33298 

0.71801 

0.12256 

0.79731 

0.24137 

0.88967 

0.26711 

0.62891 

0.35930 

0.79482 

0.10249 

0.81426 

0.28627 

0.90497 

0.38574 

0.79464 

0.44870 

0.92421 

0.30004 

0.94257 

0.39389 

1.01829 

0.23641 

1.33182 

0.08885 

1.36057 

Table  3.4.  Reconstruction  Error  of  the  CMW  Using  90%  and  95%  Compression  and  Varying  the 
Number  of  Bits  to  Represent  the  Magnitude  and  Phase 


The  number  of  phase  bits  were  decreased  for  some  experiments  because  speech  signals  are 
often  not  as  sensitive  to  the  phase  content  of  the  transform  coefficients.  This  approach  helped  to 
improve  the  signal  reconstruction  without  increasing  the  bit  rate.  For  example,  using  5  bits  for 
magnitude  and  3  bits  for  phase  (5/3)  produced  a  better  sounding  reconstruction  (and  lower  numer¬ 
ical  errors)  than  using  4  bits  for  both  the  magnitude  and  phase.  By  comparing  the  reconstructed 
signal  to  the  required  bic  rate,  the  different  coding  schemes  could  be  compared.  The  techniques 
were  compared  back  to  the  4  bits  magnitude/4  bits  phase  (4/4)  case  at  95%  compression;  signal 
reconstructions  that  were  clearly  worse  and  techniques  requiring  a  higher  data  rate  were  thrown 
out.  Two  of  the  quantization  schemes  appear  to  be  superior  to  the  other  methods  tested;  the 
5/3  scheme  at  95%  compression  and  the  3/2  technique  at  90%  compression  have  better  sounding 
reconstructed  signals  than  the  baseline  4/4  technique  and  also  require  a  lower  data  rate  to  achieve 
the  reconstruction. 

Note  that  the  numerical  error  values  do  not  always  point  toward  the  best  audio  reconstruc¬ 
tion  (25).  The  Li  values  are  not  consistent  with  the  RMS  values  and  neither  criteria  always  points 
to  the  higher  quality  speech  signal.  The  3/2  (90%)  case  has  a  much  higher  RMS  error,  but  a 
lower  Li  error  than  the  5/3  (95%)  case.  The  reconstructed  quality  and  intelligibility  of  the  3/2 
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(90%)  signal  is  slightly  better  than  the  5/3  (95%)  reconstruction.  However,  the  5/3  scheme  requires 
slightly  fewer  bits  per  second  than  the  3/2  case.  The  original  speech  signal  and  the  5/3  and  3/2 
reconstructions  are  plotted  in  figures  3.5  and  3.6. 


Position  Information:  The  spectral  position  information  can  be  coded  and  transmitted  in 
many  different  ways.  Two  methods  were  evaluated: 

1.  Sending  the  position  information  with  each  individual  non-zeroed  coefficient. 

2.  Transmitting  a  large  positional  word  prior  to  each  window  that  tells  where  the  non-zeroed 
coefficients  belong. 

The  first  method  requires  enough  bits  for  each  coefficient  to  account  for  every  possible  position.  A 
128  point  window  requires  a  seven  bit  word  for  each  coefficient  (27  =  128  possible  positions)  (2). 
The  second  method  transmits  the  same  information,  but  sends  it  prior  to  the  window  instead  of 
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0>) 

Figure  3.6.  Reconstruction  of  Sentence  Using  (a)  3  Mag  Bits/2  Phase  Bits  at  90%  Compression 
and  (b)  5  Mag  Bits/3  Phase  Bits  at  95%  Compression 
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with  each  coefficient.  A  128  point  window  would  have  a  128  bit  word  that  would  contain  a  1  at 
the  positions  where  coefficients  were  being  transmitted.  At  compression  ratios  of  90%  or  more,  the 
first  method  generally  requires  the  fewest  number  of  bits  because  of  the  low  number  of  coefficients 
which  cure  being  transmitted.  The  relationship  between  the  length  of  the  window  and  the  amount  of 
compression  demonstrates  which  method  to  use.  Method  two  always  sends  the  number  of  positional 
bits  equal  to  the  length  of  the  window,  while  the  first  method  sends  (1  —  c)L  log2  L  bits  (where  L  is 
the  window  length  and  c  is  the  percent  compression).  Therefore,  for  a  given  length  of  window,  if: 


c  >  1  — 


1 

log2  L 


(3.16) 


then  it  is  most  economical  to  transmit  the  explicit  position  with  each  term  (method  one).  If  the 
compression  is  less  than  the  value  given  in  equation  3.16,  then  the  positional  vector  should  be 
transmitted.  For  windows  of  length  128,  the  positional  vector  would  not  be  more  advantageous 
until  a  compression  ratio  of  85%  or  lower  was  used. 

The  first  method  still  requires  almost  half  of  the  system’s  bit  rate  to  be  used  just  for  positional 
information.  The  window  can  be  broken  down  into  smaller  frames  to  reduce  the  amount  of  overhead. 
For  large  compression  ratios,  there  will  be  many  zeros  within  this  information  word.  Instead  of 
transmitting  7  bits  to  give  the  exact  position  of  the  coefficient,  the  position  difference  between 
two  coefficients  could  be  transmitted.  This  difference  would  be  limited  to  a  portion  of  the  larger 
window.  For  example,  the  128  point  window  can  be  broken  into  eight  16  point  sub-windows.  If 
there  are  no  coefficients  within  16  samples,  then  a  one  bit  “dummy”  coefficient  (equal  to  zero)  is 
inserted  to  hold  the  16th  samples  place.  Thus  the  total  bits  sent  prior  to  the  window  can  be  reduced 
with  a  “divide  and  conquer”  approach  (32)  (35:4-4).  The  maximum  number  of  “dummy”  bits  to 
break  a  128  sample  window  into  16  sample  sub-windows  is  seven  (need  a  “dummy”  at  only  one  end 
point).  The  symmetry  of  the  complex- valued  Malvar  wavelet  coefficients  can  be  used  to  reduce  the 
amount  of  positional  information  required.  The  position  can  be  conveyed  with  a  six  bit  word  (64 
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sample  window  with  symmetry).  This  can  be  transmitted  with  eight  eight-point  sub-windows  and 
only  four  “dummy”  bits  will  be  required  for  the  64  sample  window. 

Examples:  The  bit  rate  for  the  transmission  of  the  spectral  coefficients  can  be  computed  by 
following  the  above  steps.  This  example  will  assume  a  sampling  rate  of  8  KHz,  128  point  windows, 
3  bits/magnitude  coefficient,  2  bits/phase  coefficient,  and  90%  compression. 

1.  8000  samples/sec  /  128  samples/window  =  62.5  windows/second. 

2.  90%  compression  of  a  128  point  window  =>  12  samples/window. 

3.  Symmetry  of  CMW  =>  6  samples/window 

4.  5  bits/sample  (magnitude  and  phase  of  coefficient) 

-1-  3  bits/sample  (coefficient  position  within  8  point  sub-window) 

8  bits/sample 

5.  4  bits/window  (“dummy”  coefficients) 

6.  (8  bits/sample)  (6  samples/window)  +  (4  bits/window)  =  52  bits/window 

7.  (62.5  windows/second)(52  bits/window)  =  3.25  Kbits/second 

The  next  example  utilizes  5  bits/magnitude  coefficient,  3  bits/phase  coefficient,  and  95% 
compression: 

1.  8000  samples/sec  /  128  samples/window  =  62.5  windows/second. 

2.  95%  compression  of  a  128  point  window  =►  6  samples/ window. 

3.  Symmetry  of  CMW  =>•  3  samples/window 

4.  8  bits/sample  (magnitude  and  phase  of  coefficient) 

+  3  bits/sample  (coefficient  position  within  8  point  sub- window) 

11  bits/sample 
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5.  4  bits/window  (“dummy”  coefficients) 

6.  (11  bits/sample)(3  samples/window)  +  (4  bits/window)  =  37  bits/window 

7.  (62.5  windows/second)(37  bits/window)  =  2.31  Kbits/second 

Even  with  the  higher  number  of  bits/coefficient,  the  higher  compression  leads  to  a  lower  data 
rate.  The  95%  compressed  version  did  not  have  the  same  quality  output  that  the  90%  version 
had,  but  the  bit  rate  is  significantly  reduced.  The  CMW  competes  very  well  with  current  speech 
encoders.  Typical  low  bit  rate  speech  encoders  used  today  such  as  the  LPC-10  coder  (used  in  the 
STU-III  secure  voice  system)  utilizes  2400  bps  (2). 

The  lapped  transform  allows  us  to  overlap  the  data  and  smooth  the  boarder  regions  without 
increasing  the  number  of  windows  transmitted  per  second.  This  is  due  to  the  fact  that  the  Malvar 
wavelet  “folds”  the  overlapped  data  back  into  the  window  before  it  calculates  the  coefficients  used  in 
the  transmission.  The  number  of  bits/second  without  compression  (but  accounting  for  symmetry) 
would  have  been: 

(8  bits/sample)(64  samples/window)(62.5  windows/second)  =  32  Kbits/sec 
The  percent  compression  of  the  overall  bit  rate  is  thus  72%  versus  the  original  95%  compression 
ratio  of  the  transform  coefficients.  This  appears  to  be  a  dramatic  decrease  in  compression,  however 
the  95%  ratio  was  simply  the  number  of  coefficients  compressed.  It  did  not  take  into  account  the 
positional  information  or  the  number  of  bits/coefficient.  The  72%  compression  achieved  still  yielded 
a  data  rate  less  than  2400  bps. 

Further  reductions  in  the  bit  rate  may  be  obtained  by  using  more  elaborate  coding  and  com¬ 
pression  techniques.  Huffman  developed  a  method  to  minimize  the  expected  number  of  transmitted 
bits  (9).  He  optimized  the  bit  pattern  by  linking  the  frequency  of  occurrence  of  a  given  value  to 
the  number  of  bits  to  code  that  value  such  that; 
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Expected  Number  of  Bits  =  (Frequency  of  Occurrence)(Bits  for  the  Values) 

The  Malvar  wavelet  has  also  been  used  to  find  voicing  probabilities  in  speech  patterns  (39).  The 
Malvar  wavelet  provides  an  automatic  segmentation  into  phonetic  units,  such  that  the  frequencies 
center  of  mass  associated  with  each  phonetic  unit  can  be  used  to  get  a  voiced/unvoiced  segmentation 
algorithm.  Voiced  and  unvoiced  speech  can  then  be  coded  and  compressed  differently  to  produce 
a  higher  quality  reconstruction  at  a  lower  bit  rate. 

Additional  Promising  Applications 

A  complex-valued  lapped  orthogonad  tramsform  would  be  a  powerful  tool  for  image  process¬ 
ing.  A  non-orthogonad  complex-valued  lapped  transform  (CLT)  has  been  defined  and  utilized  by 
Young  and  Kingsbury  to  perform  frequency-domain  motion  estimation  (41).  Young  and  Kingsbury 
extended  the  lapped  tramsform  to  form  a  complex  transform  because  data  shifts  in  the  spatial 
domain  have  more  meaningful  interpretations  in  the  complex  domaun  (primarily  as  phase  shifts  of 
the  complex  coefficients)  (41).  The  CLT  gives  smoother  motion  fields  than  the  traditional  block 
matching,  due  to  the  use  of  overlapping  windows  with  no  block  edge  discontinuities.  The  complex¬ 
valued  lapped  orthogonal  transform  proposed  by  Suter  amd  Oxley  may  yield  even  clearer  motion 
fields  because  of  its  perfect  reconstructive  properties. 

There  are  many  image  processing  applications  which  could  be  improved  through  the  use  of  a 
complex  lapped  orthogonal  tramsform,  the  CMW.  The  MITRE  Corporation  has  proposed  am  ob¬ 
jective  image  quadity  measure  that  is  highly  dependent  upon  the  reconstructive  capabilities  of  the 
short  time  Fourier  tramsform  (20).  This  quadity  measure  taikes  the  two  dimensionad  Fourier  trans¬ 
form  of  an  image,  bamdpass  filters  it,  and  reconstructs  the  image  with  amother  Fourier  transform. 
Potential  applications  in  the  image  processing  field  are  important  for  the  complex  Madvar  wavelet. 

RADAR  signad  amalysis  is  another  field  that  could  take  advantage  of  am  orthogonad  complex¬ 
valued  transform.  Phase  information  is  very  important  to  the  study  of  RADAR  signals  and  the 
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Fourier  transform  is  the  basic  tool  for  analysis  of  the  signals.  Many  spectral  estimation  techniques 
are  used  to  improve  the  resolution  and  accuracy  of  the  FFT.  The  complex- valued  Malvar  wavelet 
could  assist  these  techniques. 

Summary 

A  complex  Malvar  wavelet,  which  maintains  orthogonality,  was  introduced.  This  algorithm 
was  programmed  in  the  DOD  standard  language  of  Ada.  The  code  was  tested  and  evaluated  to 
demonstrate  that  the  transform  produces  perfectly  reconstructed  signals.  The  CMW,  with  high 
compression  ratios,  was  shown  to  be  able  to  produce  intelligible  speech  and  to  compete  with  current 
low  bit-rate  vocoders.  Imaging  and  RADAR  processing  techniques  for  the  CMW  were  introduced 
as  possible  applications. 
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IV.  Real-Valued  Malvar  Wavelets 


Introduction 

The  sinusoidal  Malvar  wavelet  (or  generalized  lapped  orthogonal  transform)  is  developed  in 
much  the  same  way  as  the  complex-valued  Malvar  wavelet,  however  the  algorithm  for  the  sinusoidal 
wavelet  is  more  complicated  because  it  is  a  real  basis  function  and  it  allows  for  variable  size  windows 
and  overlaps.  The  complete  algorithm  is  presented  in  (34)  and  (26).  The  benefit  of  using  this  real¬ 
valued  basis  is  that  the  coefficients  will  be,  by  definition,  real.  Therefore  if  the  same  amount  of 
information  can  be  sent  with  one  real  coefficient  vice  two  coefficients  (magnitude  and  phase),  there 
is  a  potential  for  a  decrease  in  the  bit  rate.  The  question  that  must  be  answered  is  whether  the 
sinusoidal  MW,  at  a  given  percent  compression,  can  produce  the  same  quality  signal  as  the  complex 
MW.  The  following  steps  summarize  the  derivation  of  the  sinusoidal  Malvar  wavelet  as  proposed 
by  Suter  and  Oxley  (34): 


1. 


Define  a  Basis.  The  basis  function  is  defined  as  the  even  and  odd  extensions  of  a  sinusoidal 
basis: 


('(*+5)  [s33y]) 

fi'k{x)  =  sin  (<k + 7)  [sjg^] ) 

\J <*y+j-7T  sin  (*(*  +  !)  [s^]) 


,  —OO  <  X  <Oj  —  Cj 

,  a j  -  Cj  <  x  <  a, 

,  Oj  <  X  <  flj  +  l 
>  aj+ 1  <  x  <  aj+ 1  +  eJ+l 
,  aj+i  +  ej+ 1  <  X  <  00. 


(4.1) 


Note  that  the  odd/even  extensions  have  been  reversed  from  the  complex-valued  case.  It  does 
not  matter  what  order  the  extensions  are  taken,  as  long  as  the  equations  remain  consistent 
throughout  the  algorithm. 

2.  Obtain  data  for  some  interval;  a,  —«;<£<  aJ+i  +  cJ+i  as  shown  if  figure  4.1. 
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dj  -  Cj  a  j  a,  +  fj 


<»j+i  -  «j+i  ai+ 1  aj+ 1  +  (j+i 


Figure  4.1.  Partitioned  Axis 


3.  Multiply  the  data  by  the  window  ( Wj(x )). 


,  -oo  <  x  <  a,-  -  (j 


8in(^-{x-[aj  -fj]»  ,  aj-(j  <x<aj  +  (j 

wj(x)=<  i  ,  aj+ej  <x<aj+i-(j  +  \ 

cos(?«7f7{*  ~  lai+1  "  ei+i3})  «  ai+1  ~  ti+1  <x^  ai+i  +  £i+i 

0  ,  Oj+1  +  £j+i  <  x  <  oo. 


4.  Fold  the  sequence;  fold  in  edges  upon  itself. 

' 

^(xjti^x)  —  s(2dj  —  x)wj(2dj  —  x)  ,  dj  <  x  <  dj  -f  (j 

Hj,l{*)  =  <  s(x)  ,  dj+Cj  <X<  dj  +  l  -  € j+1  (4-3) 

s(x)t«j(x)  +  s(2aJ+i  -  x')Wj{2dj+i  -  x)  ,  aJ+1  -  e;+i  <  x  <  dj+i 


5.  Define  a  new  array. 


Pj.i  =  sin  /°r  l  =  0,l,  -,Nj 


6.  Define  the  even  extension  of  0jj. 


Pi,l  = 


Pj,i  ,  0<l<Nj-i 
PjjNj-i  ,  Nj<l<2Nj-l 
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7.  Perform  an  FFT  of  length  2 Nj  on  this  even  extended  sequence. 

0.iei2,kt/iN}\  (4.6) 


8.  Interpret  the  results  of  the  FFT. 


Ofj,  0  —  y/2NjDjto 

(4.7) 

&j,k  —  aj,k—i  “1"  2^/2ATj  Djtt 

(4.8) 

9.  Define  the  odd  extension  of  the  spectral  coefficients. 


ai.  * 


0  <  *  <  Nj  -  1 

Nj  <  k  <  2  Nj  -  1 


(4.9) 


10.  Perform  an  Inverse  FFT  of  length  2 Nj  on  and  define. 


Hjj  =  \j2NjIm 


,  2Nj-l 

—  y  a jte^- 

k  1 


(4.10) 


11.  Reconstruct  the  signal  by  unfolding  the  data  and  dividing  by  the  window  to  fold  the  edges 
back  out. 


Sj,i  = 


^ _ ii!  x  ‘  i  —  a  • 

Wj(2aj  -  Xjil)Hj-iiNj_1-.i  +  Wj(xjjHjti)  ,  aj  <  Xjj  <  aj  +  (j 

Hj,i  i  aj  +  tj  <  xj,i  <  aj+i  —  <j+i 

wi(xi,i)Nj,i  -  Wj(2aj+i  -  *i.0^+i,jvi-i  ,  aj+i  -  tj+ i  <  xj,i  <  Oj+i 


(4.11) 
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Validation 


The  first  step  in  validating  the  proposed  system  is  to  ensure  that  the  basic  Malvar  wavelet 
code  is  correct.  The  program  written  by  Raduenz  was  not  perfect  because  he  was  concerned  with 
only  integer  values  of  signals  and  the  prime  emphasis  of  his  work  was  not  on  comparing  the  results 
of  different  systems.  The  Ada  code  was  significantly  modified  so  that  it  would  give  accurate  results 
for  the  comparison  of  multiple  systems.  The  complete  program  is  provided  in  Appendix  C.  The 
important  modifications  include: 

(1)  Allow  the  system  to  read  in  and  write  out  floating  point  values  (with  accuracy  to  the  12th 
decimal  place). 

(2)  Enable  the  system  to  perform  finite  length  processing  without  having  to  duplicate  data  or  pro¬ 
cess  intervals  multiple  times.  This  was  achieved  by  “wrapping”  the  end  points  of  the  signal,  so  that 
the  first  value  was  considered  as  the  point  immediately  after  the  last  value  of  the  data  set. 

(3)  The  starting  and  stopping  points  were  adjusted  to  correctly  match  the  algorithm  for  equa¬ 
tions  4.10  and  4.11. 

(4)  The  percent  compression  and  error  calculations  code  (discussed  in  Chapter  III  and  listed  in 
Appendix  E)  were  added  to  the  baseline  program. 

The  original  code  produced  normalized  RMS  errors  of  about  10-5  when  reproducing  a  simple 
unit  step  or  sinusoid.  The  reproduction  of  a  unit  step  is  shown  in  figure  4.2a.  This  shows  that  the 
overall  reproduction  is  very  close,  but  not  perfect  and  that  the  final  window  is  being  distorted  by 
the  algorithm.  The  modified  program’s  output  is  plotted  in  figure  4.2b.  This  is  almost  exactly  equal 
to  1.00. .0  and  has  no  distortions  in  the  overlaps.  The  normalized  RMS  error  for  the  reconstruction 
is  just  over  10“ 13.  The  reconstruction  errors  for  three  different  signals  are  given  in  table 


Input  Data 

Relative  L 2 

Loo 

RMS 

MSWOUOM 

Unit  Step 

Sinusoid 

Utterance  of  “DARK” 

m 

■  H 

jgjfB 

Table  4.1.  Reconstruction  Errors  Associated  with  the  Real- Valued  Malvar  Wavelet 


Compression 

Raduenz’s  work  suggested  that  the  homomorphic  filtering  process  may  be  superior  to  the 
use  of  the  sinusoidal  Malvar  wavelet  (LOT)  alone  (26:4-30).  The  homomorphic  filtering  allowed 
Raduenz  to  obtain  better  reconstruction  of  speech  signals  at  comparable  compression  ratios.  The 
goal  of  this  research  is  to  verify  those  results,  investigate  other  homomorphic  filters  to  determine 
if  there  is  a  best  method  for  a  given  signal,  and  then  explain  why  a  particular  method(s)  worked 
best. 

This  chapter  outlines  how  the  problem  was  attacked.  It  outlines  the  different  types  of  homo¬ 
morphic  filters  evaluated,  the  technique  used  to  threshold  (compress)  the  resulting  coefficients,  the 
error  criteria  used  to  evaluate  the  results,  and  finally  the  results  of  the  experiments. 


Homomorphic  Filters 

In  its  most  general  form,  a  homomorphic  system  can  be  represented  by  algebraically  linear 
(homomorphic)  mappings  between  the  input  and  output  signals  (22:771).  This  system  can  then 
be  decomposed  into  three  separate  homomorphic  systems;  a  characteristic  system,  a  linear  system, 
and  the  inverse  characteristic  system. 

There  were  three  distinct  systems  used  in  this  research.  Only  one  of  them  was  technically  a 
homomorphic  system.  System  1  is  shown  in  figure  4.3.  The  generalized  Malvar  wavelet  (GMW) 
acts  as  the  characteristic  system  and  five  different  functions  were  used  as  the  linear  system.  The 
five  linear  filters  (and  their  inverses)  that  were  utilized  include: 
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Figure  4.3.  System  1:  Homomorphic  System  Utilizing  the  Malvar  wavelet 

1.  No  Linear  Transform:  lx 

2.  Square  Root:  (Square) 

3.  Natural  Logarithm:  (Exponential) 

4.  10*  +  1:  ((*  -  1)/10) 

5.  Arc  Hyperbolic  Sine:  (Sinh) 

where,  hyperbolic  sine  and  arc  hyperbolic  sine  are  defined  as  (31): 


sinh(z)  =  0.5*(e*-e  *)  (412) 

sinh_1(x)  =  ln(x  +  y/x2  +  1)  (4-13) 


The  three  systems  consist  of  an  analysis  bank,  which  transforms  the  coefficients  prior  to 
transmission,  and  a  synthesis  bank,  where  the  transmitted  coefficients  are  reconstructed  back  into 
the  original  signal.  For  System  1,  the  analysis  bank  is  composed  of  the  forward  GMW  followed 
by  one  of  the  given  linear  filters  and  then  an  inverse  GMW.  The  coefficients  from  this  bank  are 
then  to  be  compressed  (thresholded  at  a  certain  level),  transmitted,  and  sent  to  the  synthesis  bank 
for  signal  reconstruction.  The  synthesis  bank  is  made  up  of  a  forward  GMW,  the  inverse  of  the 
linear  filter  used  in  the  analysis  bank,  and  finally  an  inverse  GMW.  System  2,  figure  4.4,  represents 
the  system  used  by  Raduenz.  This  system  is  comprised  of  the  GMW,  a  linear  filter,  and  another 
forward  GMW.  The  synthesis  filter  is  constructed  from  the  inverse  GMW,  an  inverse  linear  filter, 
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Figure  4.4.  System  2:  Homomorphic- Like  System  Utilising  the  GMW 


Figure  4.5.  System  3:  Malvar  Wavelet  with  Linear  Filter 

and  finally  another  inverse  GMW.  The  final  system  (system  3)  used  in  these  experiments  consists 
only  of  the  GMW  followed  by  a  linear  filter  as  shown  in  figure  4.5.  There  are  15  different  systems 
to  be  analyzed;  five  different  linear  filters  in  each  system.  Raduenz  performed  the  homomorphic 
filtering  by  completing  the  first  forward  GMW,  taking  the  natural  logarithm  of  all  the  GMW 
coefficients,  and  then  performing  another  forward  GMW.  In  order  to  be  able  to  perform  the  linear 
filtering  in  real-time,  a  procedure  was  added  to  the  GMW  program  to  implement  the  linear  filter  on 
each  window  of  data  after  the  characteristic  system  (GMW)  had  transformed  that  given  window. 
The  linear  filter  will  then  pass  the  set  of  data  on  to  the  next  function  (Inv  GMW  for  System  1). 
To  accomplish  this,  a  step  was  added  to  the  procedure  Do  Work  of  the  Lapped  Transform  program 
that  calls  the  new  procedure  Linear  Transform.  The  steps  which  comprise  the  Linear  Transform 
procedure  in  the  Lapped  Transform  program  are  found  in  Appendix  D. 

The  first  time  the  Linear-Transform  procedure  is  called,  the  value  of  Inverse  is  set  to  zero  so 
that  the  forward  transform  is  performed.  Subsequent  times  the  procedure  is  called,  Inverse  is  set 
to  one,  so  that  the  inverse  transform  is  used.  The  type  of  linear  transform  (log,  y/z,  etc)  is  selected 
by  the  user,  prior  to  the  start  of  any  processing.  Note  that  the  phase  of  the  GMW  coefficients  do 


4-8 


not  have  to  be  tracked  because  the  coefficients  of  the  GMW  (as  defined  at  the  beginning  of  this 
chapter)  are  always  real- valued.  However  the  positive  and  negative  values  could  be  expressed  as  a 
+/—  or  0°/180°  phase.  The  negative  values  pose  a  problem  for  the  logarithm  and  the  square  root 
functions.  This  could  be  corrected  for  by  expressing  the  functions  with  a  phase  term  and  utilizing 
the  complex  function,  log[X(f)]  =  Zn|X(/)|  +  j/.X(f).  In  order  to  avoid  any  kind  of  phase  term, 
the  coefficients  can  be  offset  so  that  the  minimum  coefficient  value  is  1.0  and  the  minimum  log 
amplitude  is  therefore  0.0.  The  offset  can  then  be  removed  after  the  inverse  natural  log  (er ).  For 
example,  if  there  exists  a  set  of  coefficients  p(x),  an  offset  of  y  —  minimum[p(z)]  +  1  is  used  with 
the  logarithm: 

1.  Coefficient  value  p(x) 

2.  Add  an  offset  of  y:  (p(x)  +  y) 

3.  Take  the  natural  log:  ln(p(x)  +  y) 

4.  Perform  additional  linear  transforms:  (Lapped  Trans  or  Inv  LOT) 

5.  Take  inverse  natural  log:  el-p<-T>,+y)  =  p(x)  +  y 

6.  Subtract  the  offset:  p(x)  +  y-  y  =  p(x)  +  (min  +  1)  —  (min  +  1)  =  p(x) 

This  same  procedure  can  be  applied  to  the  square  root,  with  the  offset  equal  to  the  minimum 
coefficient  value,  because  the  square  root  of  zero  is  defined.  This  method  works  well  when  utilizing 
all  of  the  coefficients.  However,  compressing  the  coefficients  (setting  coefficients  to  zero)  introduces 
a  non-linearity  to  the  system.  This  non-linearity  means  that  the  minimum  valued  being  carried 
through  the  system  will  distort  the  reconstructed  signal.  A  third  method  can  be  used  to  track  the 
sign  of  the  original  coefficients.  The  sign  (+/-)  of  the  original  coefficient  is  stored  in  a  vector  and 
then  the  absolute  value  of  the  coefficient  is  used  with  the  function.  The  sign  is  then  put  back  on 
the  signal  after  reconstruction  (e*  or  x2).  The  benefit  of  this  method  is  that  the  sign  vector  will 
not  distort  the  resulting  signal  when  coefficients  are  eliminated. 
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System  Validation 

Each  complete  system  must  now  be  validated  to  ensure  packages  are  working  together  prop¬ 
erly.  To  validate  each  program  and  linear  transform  system,  known  signals  were  passed  through  the 
systems  with  0%  compression.  Perfect  reconstruction  is  expected  when  there  is  no  lose  of  informa¬ 
tion  (no  coefficients  thresholded).  Allowing  for  some  round-off  error,  each  filter  and  linear  system 
produced  perfectly  reconstructed  signals.  The  speech  signals  sounded  identical  and  reconstruction 
errors  were  very  low.  The  numerical  error  criteria  discussed  in  Chapter  III  were  used  to  compare 
the  results. 


Numerical  Error  Tests:  The  reconstruction  errors  (£21  £».  RMS,  and  RIV)  were  collected 
and  compared.  System  1  (GMW,  Linear  Function,  Inv  GMW:  see  figure  4.3)  and  system  2  (GMW, 
linear  Function,  GMW:  see  figure  4.4)  produced  nearly  identical  output  errors.  There  are  many 
different  ways  in  which  both  systems  may  be  coded  as  part  of  the  lapped  transform  routine,  however 
they  always  produced  very  similar  levels  of  reconstruction  error.  System  3  (GMW,  Linear  Function: 
see  figure  4.5)  also  provided  near  perfect  reconstruction  at  0%  compression  (using  125  point  windows 
with  50  point  overlap).  The  numerical  errors  have  been  tabulated  below  for  a  2000  point  sinusoidal 
input. 


Linear  Function 

Relative  L3 

Loo 

RMS 

RIV  (1.0  +  ) 

1 

2.284x10"* 

6.037xl0"lu 

5.772xl0”12 

-7.555x10"“ 

v? 

3.058x10-* 

7.981xlO-10 

8.626xl0”12 

-1.354x10- 10 

log(x) 

2.327xl0“5 

5.052xl0"10 

5.283xl0-12 

—7.840x10” 11 

10x+  1 

2.284x10"* 

6.090xl0"10 

5.818xl0-13 

-7.555x10- 11 

sinh(x) 

2.915x10-* 

7.261xlO"10 

7.736x10" 12 

-1.230x10-“ 

Table  4.2.  Reconstruction  Errors  Associated  with  System  1 


All  four  functions  (lx,  y/x,  log(x),  10z+l,  and  sinh(x))  produced  identical  outputs  for  system 
3.  These  reconstructed  signals  had  lost  very  little  information  and  were  perfect  reconstructions  of 
the  originals.  The  output  signals  produced  by  system  1  and  2  were  not  as  well  defined.  The  linear 
functions  (lx,  and  lOx  +  1)  produced  nearly  identical  output  signals.  Note  that  these  functions 


4-10 


Linear  Function 

Relative  L? 

loo 

RMS 

RIV  (1.0  +  ) 

1 

2.464zl0_s 

6.525zl0-10 

5.876zl0~12 

-8.795Z10"11 

V? 

3.081Z10"5 

6.472zlO~10 

7.778zl0-12 

—  1.374zl0_1° 

log(z) 

5.932zl0-5 

1.954zl0-9 

2.800zl0-11 

— 5.094zl0_1° 

lOz  +  1 

2.464zl0-5 

6.550zl0-10 

5.880zl0-12 

— 8.795zl0-11 

sinh(z) 

3.396zl0"5 

7.532zl0_1° 

9.419zl0"12 

— 1.670zl0_1° 

Table  4.3.  Reconstruction  Errors  Associated  with  System  2 


Linear  Function 

Relative  L? 

loo 

RMS 

RIV  (1.0  +  ) 

l 

1.615zl0_a 

2.886zl0“ ia 

-3.777zl0-“ 

1.615zl0-5 

EC  2  29 

2.886zl0-12 

-3.777zl0~n 

log(z) 

1.615zl0-5 

2.886zl0-12 

— 3.777zlO“n 

1.615zl0"5 

2.886zl0-12 

— 3.777zl0"11 

sinh(z) 

1.615zl0-5 

ebbssh 

2.886zl0-12 

— 3.777zl0-11 

Table  4.4.  Reconstruction  Errors  Associated  with  System  3 


also  had  the  lowest  reconstruction  errors  for  system  1  and  2.  The  other  functions  produced  slightly 
higher  errors  than  the  lz  and  lOz  +  1  functions.  The  recorded  errors  are  so  small  that  they  are 
very  sensitive  to  the  number  of  flops  (multiplies  and  adds)  performed  by  the  system.  System  3,  the 
most  basic  procedure,  uses  the  least  number  of  flops  and  thus  has  the  lowest  reconstruction  errors. 
System  2  has  the  greatest  number  of  flops  because  of  the  way  the  system  had  to  be  constructed.  The 
original  lapped  transform  program  was  not  written  to  have  two  forward  transforms  in  succession, 
so  it  has  more  coefficient  calculations  than  do  the  other  two  systems.  It  follows  then  that  system 
2  has  slightly  higher  reconstruction  errors. 

All  of  the  reconstruction  errors  are  extremely  small  (the  largest  normalized  RMS  error  = 
2.800zl0-u)  and  easily  demonstrate  the  perfect  reconstruction  properties  of  the  systems. 

Graphical  and  Audio  Tests:  Graphical  and  audio  tests  were  also  made  to  measure  the  signal 
reconstructions.  These  again  showed  how  well  all  of  the  systems  were  restored  when  there  is 
no  thresholding  of  the  coefficients.  Figure  4.6  shows  the  sinusoid  that  was  used  for  the  initial 
testing  of  the  systems.  Figure  4.7  (a)  is  the  output  of  system  3  using  no  linear  function  (GMW, 
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Inverse  GMW).  The  output  of  the  systems  using  other  functions  (log,  sinh,  etc)  were  identical  to 
figure  4.7.  Figure  4.7  (b)  shows  the  output  from  system  1  using  no  linear  function  (GMW,  Inverse 
GMW,  GMW,  Inverse  GMW).  This  restoration  is  also  perfect.  The  outputs  of  the  other  functions 
performed  virtually  the  same  as  the  lx  case. 


Figure  4.6.  Original  Sinusoid  Signal 


After  testing  the  system  with  a  step  function  and  a  sinusoid,  speech  files  from  the  TIMIT 
database  were  used.  Words  and  parts  of  words  that  were  2000  points  long  were  tested  with  systems 
1,  2,  and  3.  This  allowed  for  more  complicated  signals  to  be  transformed  by  the  programs  and  to 
actually  listen  to  the  results  instead  of  simply  looking  at  numbers  and  plots.  Figure  4.8  shows  the 
graphical  display  of  the  utterance  for  the  word  “dark.”  The  reconstructions  utilizing  system  3,  and 
system  1  are  plotted  in  figures  4.9(a)  and  (b). 

The  three  plots  of  the  word  “dark”  look  very  much  the  same,  but  that  does  not  guarantee 
that  they  will  sound  identical.  When  played  through  the  Ariel  audio  interface,  the  words  could 
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Figure  4.7. 


M 


(b) 

(a)  GMW  (System  3)  and  (b)  Homomorphic  (System  1)  Reconstruction  of  the 
Sinusoid 
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Figure  4.8.  Utterance  of  the  Word  "Dark” 

not  be  differentiated.  However,  short  one  syllable  words  often  sound  alike  because  there  is  so  little 
time  for  the  ear  to  differentiate  them.  The  next  step  was  to  utilize  entire  sentences  of  speech.  The 
results  from  transforming  entire  speech  files  was  very  similar  to  the  previous  results.  All  of  the 
systems  provided  perfect  reconstruction  using  all  criteria;  numerical  values,  plotting,  and  listening. 
Listening  provided  the  least  differentiation  between  the  systems.  The  utterance,  “She  had  your 
dark  suit  in  greasy  wash  water  all  year,”  was  reconstructed  using  the  three  systems.  All  provided 
perfect  reconstruction  of  the  speech  signed  (the  worst  normalized  RMS  error  was  10~9).  All  of  the 
signeds  sounded  identiced  when  played  via  xwaves.  There  were  no  added  noise,  pops,  or  spikes  in 
the  recreated  signals. 

Observations;  The  previous  section  demonstrates  that  perfect  reconstruction  was  possible 
using  homomorphic  filtering  techniques  with  Malvar  wavelets.  The  results  also  demonstrates  that 
the  most  basic  system,  consisting  of  only  the  GMW  and  the  Inverse  GMW,  yields  the  best  results 


samples 


Figure  4.11.  (a)  System  2  (log)  (b)  System  3  (sinh)  Reconstruction  of  the  Sentence  “She  had  your 

dark  suit  in  greasy  wash  water  all  year” 
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at  zero  compression.  This  is  due  to  the  fact  that  the  basic  system  had  the  fewest  multiplies 
and  additions.  Another  system  may  show  much  greater  promise  when  there  is  a  high  degree  of 
compression  of  the  coefficients. 

To  observe  the  effects  of  compression,  error  criteria  was  collected  using  10  different  data 
sets.  The  data  consisted  of  five  different  sinusoid  signals  and  five  different  speech  signals.  The 
sinusoids  varied  in  frequency  and  complexity,  while  the  speech  segments  were  short  words  and 
sounds  consisting  of  both  voiced  and  unvoiced  speech:  “dark,”  “she,”  “wash,”  “s,”  and  “ear.” 
All  ten  signals  are  plotted  in  Appendix  A.  The  amount  of  compression  was  also  varied  for  each 
signal:  0%,  80%,  and  90%  compression  were  used  for  each  data  set.  This  was  done  to  ensure  that 
a  particular  system  does  well  at  both  high  compression  rates  as  well  as  at  zero  threshold.  A  total 
of  30  data  files  were  used  for  each  transform  system  (each  of  the  ten  data  sets  compressed  to  three 
different  levels).  A  number  of  variables  were  held  constant  throughout  the  experiments,  including: 
the  data  length  (1920  points),  window  size  (128  points),  the  amount  of  overlap  (1  point),  and  the 
basis-set  (sinusoidal  basis  described  in  Chapter  II). 

Tables  4.5,  4.6,  and  4.7  list  the  average  reconstruction  errors  for  systems  1,  2,  and  3  respec¬ 
tively.  The  variance  of  each  of  the  error  criteria  are  recorded  in  tables  4.8,  4.9,  and  4.10. 


Linear  Function 

Relative 

Lqq 

RMS 

RIV 

1 

0.375381 

161.326 

1.27688 

0.942922 

lOz  +  1 

0.398733 

182.457 

1.10785 

0.939677 

v/£ 

0.413016 

192.152 

1.03613 

0.933322 

Iog(z) 

0.374349 

428.568 

1.92029 

0.964840 

sinh(z) 

3.22121 

7757.98 

31.1507 

0.993657 

Table  4.5.  Average  Errors  Associated  with  System  1 


The  results  of  the  450  experiments  show  that  the  system  which  consistently  performed  the 
best  was  system  3  (with  the  non- logarithmic  functions;  1,  10a:  +  1,  and  y/x).  At  every  level  of 
compression,  with  every  type  of  signal,  the  non-Iogarithmic  system  3  reconstructed  the  signals  with 
the  lowest  numerical  errors.  It  is  expected  that  these  functions  would  be  identical  for  system  3, 
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Linear  Function 

Relative  Li 

Loo 

RMS 

RIV 

1 

0.363028 

213.032 

1.24357 

0.950170 

lOz  +  1 

0.398701 

201.787 

1.09956 

0.948737 

0.406773 

178.925 

1.01884 

0.940515 

log(z) 

0.400874 

370.767 

1.74659 

0.956094 

sinh(z) 

2.97311 

6699.29 

27.8690 

0.996466 

Table  4.6.  Average  Errors  Associated  with  System  2 


Linear  Function 

Relative  Li 

Loo 

RMS 

RIV 

1 

0.104472 

213.757 

0.685003 

0.998475 

lOz  +  1 

0.104472 

213.757 

0.685003 

0.998475 

v? 

0.104472 

213.757 

0.685003 

0.998475 

log(z) 

1.03833 

215.032 

0.687536 

1.14262 

sinh(z) 

1.03833 

215.032 

0.687536 

1.14262 

Table  4.7.  Average  Errors  Associated  with  System  3 


because  the  system  (GMW,  L  — *  Inv  L,  Inv  GMW)  performs  the  inverse  functions  immediately 
after  the  forward  function  (with  only  compression  in  between).  The  logarithms  (ln(z)  and  sinh(z)) 
have  a  higher  distortion  across  the  compressions.  The  reason  for  this  is  that  the  natural  logarithms, 
due  to  the  nature  of  the  function,  are  thresholding  some  higher  coefficients  while  keeping  smaller 
ones  (for  example  |log(99)|  <  |log(.01)|).  This  distortion  could  be  avoided  by  adding  a  term  to 
the  coefficients  to  ensure  a  minimum  value  of  1.0  but,  as  stated  previously,  this  introduces  another 
distortion  when  the  coefficients  are  compressed.  The  sinh  function  works  relatively  well  (same 
as  the  In  function)  with  system  1,  but  it  distorts  the  signal  quite  badly  when  there  is  a  second 
transform  following  the  sinh  (systems  1  and  2). 


Linear  Function 

Relative  Li 

Loo 

RMS 

RIV 

1 

0.123228 

81316.5 

4.56385 

8.22438zl0~3 

lOz  4-  1 

0.125366 

98737.7 

4.36783 

5.35678zl0~3 

v/? 

0.128805 

114023. 

3.16373 

1.15164zl0~2 

log(z) 

0.104304 

584534. 

11.4624 

6.07249zl0~3 

sinh(z) 

21.1635 

2.37950zl0+8 

4027.31 

4.16942zl0-2 

Table  4.8.  Variance  of  Errors  Associated  with  System  1 
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Linear  Function 

e 

RMS 

RIV 

1 

0.114239 

132930. 

4.28904 

6.35135zl0-3 

lOx  +  1 

0.116201 

122780. 

4.01349 

7.76930zl0-3 

s/i 

0.122680 

102405. 

3.11571 

8.78877xl0-3 

log(z) 

0.124123 

438373. 

9.30463 

1.21154xl0-2 

sinh(x) 

18.6011 

1.94071xl0+s 

3556.47 

3.73096xl0-2 

Table  4.9.  Variance  of  Errors  Associated  with  System  2 


Linear  Function 

Relative  L2 

Loo 

RMS 

RIV 

1 

0.0125197 

151753. 

1.44600 

3.89581x10"® 

lOx  +  1 

0.0125197 

151753. 

1.44600 

3.89581x10-® 

V? 

0.0125197 

151753. 

1.44600 

3.89581x10"® 

log(x) 

5.73403 

148896. 

1.42399 

0.162145 

sinh(x) 

5.73403 

148896. 

1.42399 

0.162145 

Table  4.10.  Variance  of  Errors  Associated  with  System  3 


System  3  outperforms  the  other  systems  in  the  L2,  RMS,  and  RIV  criteria.  The  value 
(the  maximum  point-wise  difference  between  the  original  and  reconstructed  signal),  alone,  is  not  a 
good  indicator  of  reconstruction  quality.  The  L<x,  value  does  indicate  whether  the  output  signal  has 
many  large  spikes  in  the  reconstruction.  The  RIV  criteria  produces  the  most  accurate  reproduction 
when  its  value  is  equal  to  one  (information  in  reproduction  =  information  in  original).  The  RIV 
value  should  generally  not  be  greater  then  one,  because  this  would  indicate  information  being  added 
to  the  signal,  however  many  of  the  systems  distorted  the  signal  badly  enough  to  cause  the  RIV 
criteria  to  be  greater  than  one.  Therefore,  the  average  RIV  can  not  be  used  to  accurately  determine 
the  quality  of  the  reconstructed  signal. 

The  variance  (see  table  4.10)  of  system  3’s  error  values  are  also  the  lowest  for  the  i2,  RMS,  and 
RIV  criteria.  The  L «,  case  is  the  one  criteria  where  system  1  has  consistently  lower  values.  System  3 
(logarithm)  produces  incongruities  and  spikes  at  the  window  boundaries  at  high  compression  ratios. 
System  1  smoothes  out  these  spikes  and  thus  reduces  the  maximum  point  wise  error.  Figure  4.12(a) 
shows  a  portion  of  the  reconstructed  “s”  using  system  3,  logarithm,  and  80%  compression.  The 
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same  data  is  reconstructed  with  system  1  in  figure  4.12(b).  The  large  spikes  in  the  first  figure  have 
been  greatly  reduced  by  the  additional  Malvar  wavelet  manipulations  of  system  1 . 

Not  only  sure  the  numerical  errors  lower  for  the  system  3  case  but  the  graphical  and  audio 
reconstructions  are  also  superior  to  the  other  methods.  The  measure  does  not  have  much 
bearing  on  the  total  sound  or  graphical  quality  of  the  reconstruction.  Some  of  the  large  spikes 
manifest  themselves  as  pops  in  the  reproduction,  but  the  overall  quality  of  the  signal  is  superior 
with  system  3. 

As  the  amount  of  compression  increases,  the  reconstruction  errors  at  the  window  boundaries 
lead  to  higher  total  reconstruction  errors.  This  can  be  seen  by  taking  a  closer  look  at  a  given 
boundary  point  for  a  signal  with  increasing  levels  of  data  compression. 

Figures  4.13  and  4.14  show  the  reconstruction  of  a  sinusoid  with  compression  levels  of  0%, 
80%,  and  98%  respectively.  Each  window  is  128  points  long,  so  the  first  window  edge  falls  at  point 
128.  The  degradation  of  the  quality  of  reconstruction  is  easily  seen  as  the  amount  of  compression 
increases.  This  decay  in  the  quality  of  reconstruction  is  possibly  due  to  the  minimal  amount  of 
overlap  (1  point)  being  used  for  these  experiments.  One  of  the  primary  benefits  of  the  Lapped 
Transform  is  being  able  to  increase  the  amount  of  overlap  without  increasing  the  amount  of  data 
that  must  be  transmitted.  This  problem  raises  two  new  questions. 

1.  Will  an  increase  in  the  amount  of  overlap  produce  a  more  accurate  recreation  of  the  original 

signal? 

2.  Is  there  an  optimal  amount  of  overlap  which  should  be  utilized? 

These  questions  will  be  addressed  and  tested  in  Chapter  V. 

These  systems  were  used  to  transform  entire  speech  files  from  the  TIMIT  database.  Plotting 
and  listening  to  the  reconstructed  speech  led  to  the  same  conclusion  as  above.  System  3  (with 
a  non-logarithmic  function)  produces  superior  reconstructed  signals  when  compared  to  the  other 
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Figure  4.12. 
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(a)  System  3  and  (b)  System  1  Reconstruction  of  the  “S”  Sound  with  80% 
Compression. 
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Figure  4.13.  Reconstruction  of  a  Sinusoid,  0%  Compression,  Window  Edge  =  Point  128 

systems  in  the  test.  The  transformed  speech  from  system  3  was  consistently  more  intelligible  and 
had  higher  reconstruction  quality. 

It  is  difficult  to  determine  why  Raduenz  found  that  a  transform  like  system  2  led  to  greater 
compression  capabilities.  He  was  not  using  a  program  that  provided  perfect  reconstruction,  so  the 
more  complicated  system  may  have  made  up  for  the  imperfections  of  the  basic  program.  He  also 
may  have  changed  to  16  KHz  signals,  which  do  give  higher  fidelity  reconstructions.  Raduenz  had 
assumed  that  changing  from  16  KHz  to  8  KHz  would  not  effect  his  results.  This  is  essentially  true 
when  there  is  no  compression  involved.  However,  at  high  compression  ratios,  the  change  from  16 
KHz  to  8  KHz  does  make  a  difference.  For  the  experiments  from  this  thesis,  all  of  the  signals 
(including  the  original  speech  files)  were  analyzed  at  8  KHz. 
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Figure  4.14.  (a)  80%  Compression  and  (b)  98%  Compression  Reconstruction  of  Sine,  Window 

Edge  =  Point  128 
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Complex  and  Real-Valued  Malvar  Wavelets 


The  Complex  and  the  Sinusoidal  Malvar  wavelets  (CMW  and  SMW)  can  now  be  directly 
compared,  since  it  has  been  determined  that  the  basic  Malvar  wavelet  produces  reconstructed 
signals  as  good  as  or  better  than  any  other  system  that  was  considered.  It  has  been  shown  previously 
that  both  of  the  wavelet  types  produce  very  accurate  signals  even  at  quite  high  compression  ratios. 
How  do  the  two  systems  compare  when  considering  reconstruction  quality  of  speech  signals  at  high 
compression  ratios  (90%  -  99%)? 

The  two  Malvar  wavelets  were  used  to  reproduce  a  speech  file  at  compression  ratios  of  90%, 
95%,  and  99%.  The  original  signal  is  shown  in  figure  4.15a.  Figures  4.16a  and  4.17a  show  how 
the  two  transforms  reconstructed  the  signal  at  90%  compression.  There  does  not  appear  to  be  too 
much  difference  except  that  the  Complex  MW  plot  reaches  fewer  of  the  peaks  of  the  signal.  This 
can  be  shown  more  easily  by  plotting  a  smaller  segment  of  the  file.  Figure  4.15b  shows  the  original 
speech  signal  from  sample  16000  to  sample  17000.  The  transform  reconstructions  of  this  segment 
are  shown  in  figures  4.16b  and  4.17b.  It  is  clear  that  both  transforms  are  loosing  information,  but 
the  complex  wavelet  has  lost  much  more  information  at  the  same  compression  ratio.  The  difference 
in  reconstruction  is  also  shown  by  the  numerical  reconstruction  errors.  Table  4.11  lists  the  Li  and 
RMS  errors  for  the  two  transforms  at  90%,  95%,  and  99%  compression.  The  sound  of  the  SMW 
speech,  at  90%  compression,  was  very  intelligible.  Much  of  the  quality  had  been  lost.  There  was 
a  mechanical  tone  of  voice  with  ringing  bell-like  artifacts.  The  speech  may  have  warbled  and  had 
a  “spacey”  sound  to  it,  but  it  was  quit  understandable.  The  CMW  signal  was  also  intelligible, 
however  it  had  added  a  constant  low-level  noise.  There  was  a  constant  low  popping  all  the  way 
through  the  file.  At  95%  compression  (plotted  in  figure  4.18),  the  complex  transform  had  many 
sudden  starts  and  stops  at  the  beginning  and  ending  of  words.  The  low-level  popping  or  clicking 
was  also  more  evident  at  95%.  The  SMW  reconstruction  contained  fewer  sudden  starts/stops  and 
did  not  have  the  same  noise  as  the  CMW  version.  More  intelligibility  had  been  lost,  at  95%,  but 
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most  of  the  words  can  still  be  picked  out.  The  spectral  coefficients  of  the  wavelets  (real  coefficients 
only)  are  shown  in  figure  4.19.  The  spectral  coefficients  follow  the  same  general  pattern,  however 
the  SMW  coefficients  are  much  higher  magnitude  and  have  a  few  more  peaks  within  the  groups 
of  coefficients.  At  99%  compression  both  systems  have  lost  virtually  all  intelligibility.  Speech 
waveforms  are  definitely  still  present  (see  figure  4.20  a  and  b),  however  they  have  been  cut  down 
and  smoothed  over  so  much  that  real  vocalization  is  not  present. 


SMW,  90% 

|  SMW,  95 % 

CMW,  95% 

SMW,  99% 

CMW,  99% 

Li 

RMS 

0.16433 

0.35393 

0.33098 

0.71379 

0.43108 

0.93077 

0.53188 

1.14551 

0.67070 

1.44693 

Table  4.11.  L 2  and  RMS  Errors  for  the  Sinusoidal  and  Complex  Malvar  Wavelets  at  90%,  95%, 
and  99%  Compression. 
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(b) 


Figure  4.16.  Sinusoidal  MW  Reconstruction  at  90%  Compression  (a)  Complete  Signal  and  (b) 
Samples  16000  to  17000. 


4-28 


Figure  4.17. 


(b) 


Complex  MW  Reconstruction  at  90%  Compression  (a)  Complete  Signal  and  (b) 
Between  Samples  16000  and  17000. 
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Figure  4.18.  (a)  Sinusoidal  and  (b)  Complex  MW  Reconstruction  at  95%  Compression. 
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Figure  4.19.  (a)  Spectral  Coefficients  from  the  SMW  and  (b)  Real  Spectral  Coefficients  from  the 

CMW. 
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Figure  4.20.  (a)  Sinusoidal  and  (b)  Complex  MW  Reconstruction  at  99%  Compression. 


Compression/Data  Rates 

The  data  rate  calculations  from  Chapter  III  can  easily  be  repeated  for  the  sinusoidal  Malvar 
wavelet.  The  data  rate  calculations  will  use  a  window  size  of  128  samples,  sample  rate  of  8  KHz, 
and  95%  compression.  The  real-valued  wavelet  produces  128  real  coefficients  for  each  transformed 
window.  This  transform,  like  the  Discrete  Cosine  Transform,  has  no  coefficient  symmetry.  The 
amplitude  of  these  coefficients  can  be  coded  using  4  bits/coefficient.  The  bits/sec  for  the  spectral 
values  are  computed  by: 

1.  8000  (samples/sec)  /  128  (samples/window)  =  62.5  windows/second. 

2.  95%  compression  of  a  128  point  window  =*■  6  samples/window. 

3.  4  bits/sample  (amplitude  for  the  real  coefficient) 

+  7  bits/sample  (coefficient  position) 

11  bits/sample 

4.  (11  bits/sample)  (6  samples/window)  =  66  bits/window 

5.  (62.5  window8/second)-(66  bits/window)  =  4.125  Kbits/second 

However,  over  60%  of  the  system’s  4125  bits  are  used  just  for  positional  information.  The 
position  of  the  coefficients  is  being  kept  by  one  128  point  word.  This  one  large  word  can  be 
broken  down  into  smaller  words.  For  large  compression  ratios,  there  will  be  many  zeros  within  this 
information  word.  Instead  of  transmitting  7  bits  (27  =  128  possible  positions)  to  give  the  exact 
position  of  the  coefficient,  the  position  difference  between  two  coefficients  could  be  transmitted. 
This  difference  could  be  limited  to  a  portion  of  the  larger  window.  For  example,  the  128  point 
window  can  be  broken  into  eight  16  point  sub-windows.  If  there  are  no  coefficients  within  16 
samples,  then  a  one  bit  “dummy”  coefficient  (equal  to  zero)  is  inserted  to  hold  the  16th  samples 
place.  Thus  the  total  bits  sent  prior  to  the  window  can  be  reduced  with  a  “divide  and  conquer” 
approach  (32)  (35:4-4).  The  maximum  number  of  “dummy”  bits  to  break  a  128  sample  window 


4-33 


into  16  sample  sub- windows  is  seven  (need  a  “dummy”  at  only  one  end  point).  The  bit  rate  is  then 
recalculated  as: 

1.  62.5  windows/second 

2.  6  samples/window 

3.  4  bits/sample  (amplitude  for  the  real  coefficient) 

+  4  bits/sample  (coefficient  position  within  sub-window) 

8  bits/sample 

4.  7  bits/window  (“dummy”  coefficients) 

5.  (8  bits/sample)(6  samples/ window)  +  (7  bits/window)  =  55  bits/window 

6.  (62.5  windows/second)(55  bits/window)  =  3.437  Kbps 

The  CMW  maintains  a  lower  data  rate  for  the  same  compression  ratio,  however  the  intelli¬ 
gibility  of  the  reconstructed  signal  is  slightly  worse.  Note  that  the  CMW  reconstructions  for  this 
section  were  done  using  the  real  and  imaginary  coefficients.  It  was  shown  previously  that  the  mag¬ 
nitude  and  phase  gave  a  better  quality  reconstruction.  Therefore  the  CMW  and  the  SMW  do  yield 
about  the  same  quality  reconstructions  at  comparable  compression  ratios  and  the  two  transforms 
are  virtually  identical  when  both  intelligibility  and  data  rate  are  taken  into  account. 

Summary 

The  experiments  performed  and  the  data  collected  demonstrate  that  the  most  accurate 
method  of  compressing  and  reconstructing  digitized  data  files  is  also  the  most  basic  method;  the 
Sinusoid  Malvar  wavelet.  This  conclusion  makes  intuitive  sense;  by  performing  just  the  forward 
GMW  and  then  compressing  the  data,  there  are  fewer  calculations  to  be  made  and  therefore  a 
higher  degree  of  accuracy  may  be  maintained.  Experiments  also  demonstrated  that  intelligible 
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speech  signals  can  be  obtained  for  3  Kbps  or  less  using  both  the  real-valued  and  complex- valued 
Malvar  wavelet. 
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V.  Overlap  Optimization  For  Real-  Valued  Malvar  Wavelets 


Introduction 

It  was  shown  in  the  previous  chapter  that  at  high  compression  rates  it  may  be  necessary  to 
increase  the  amount  of  overlap  between  adjacent  intervals  in  order  to  smooth-out  the  transition. 
There  is  no  benefit  to  increasing  the  amount  of  overlap  for  small  amounts  of  compression,  because 
the  reconstruction  is  nearly  flawless  as  it  is.  However,  when  the  signal  is  being  compressed  80%  to 
99%  ,  increasing  the  overlap  can  improve  the  reconstruction  greatly.  The  advantage  of  increasing 
the  amount  of  overlap  is  shown  in  figures  5.1(a)  and  (b). 

Figure  5.1(a)  is  the  reconstruction  of  a  sine  wave  after  98%  compression  of  the  transform 
coefficients.  The  boundary  between  blocks  has  been  severely  distorted.  The  next  plot  shows  the 
same  signal  with  the  same  amount  of  compression,  but  with  50%  (64  point)  overlap.  The  distortion 
has  been  eliminated.  The  question  that  needs  to  be  answered  is  whether  it  is  necessary  to  always 
maximize  the  amount  of  overlap  or  is  there  some  threshold  where  continuing  to  increase  the  overlap 
yields  limited  returns. 

Selecting  An  Error  Criteria 

Some  sort  of  numerical  error  criteria  must  be  used  in  order  to  be  able  to  have  a  statistically 
meaningful  confidence  interval  for  the  amount  of  overlap.  Chapters  III  and  IV  showed  that,  in 
general,  the  lower  the  error  value  the  better  the  reconstruction.  Table  5.1  shows  how  the  error 
criteria  of  Chapters  III  and  IV  vary  with  increasing  compression  for  a  given  signal. 


Percent  Compression 

RIV 

Relative  L i 

RMS 

Squared 

0% 

0.99997 

0.01263 

6.592arl0-Y 

0.99984 

99% 

0.94209 

0.52631 

4.168xl0-3 

0.72582 

— oo 

1.00000 

7.954*10"3 

0.00000 

Table  5.1.  Enror  Criteria  Comparison 
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Figure  5.1. 


(a) 


(b) 

(a)  1  Point  Overlap  and  (b)  50%  Overlap  Reconstructions  of  a  Sinusoid,  98%  Coeffi¬ 
cient  Compression 
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Results  from  the  previous  chapters  demonstrate  that  the  RIV  and  the  criteria  are  not 
very  good  reconstruction  predictors.  The  Renyi  Information  Value  can  be  skewed  because  of  the 
non-linearity  of  the  compression.  The  RIV  is  also  always  very  close  to  1.0  and  therefore  does  not 
differentiate  between  the  reconstructions  very  well.  The  Leo  value  does  give  information  on  the  size 
and  magnitude  of  extraneous  spikes  in  the  reproduced  signal.  These  spikes  are  not  always  to  worst 
part  of  the  reproduction.  Some  signals  were  reproduced  with  large  spikes,  but  were  otherwise  very 
close  to  the  original  signed;  while  other  reconstructions  were  very  smooth  signals  but  did  not  match 
the  original  very  well. 

The  two  error  criteria  that  do  follow  the  reconstruction  fairly  well  are  the  relative  Li  and  the 
RMS  errors.  The  Li  was  selected  because  it  is  normalized  between  zero  and  one  (see  table  5.1). 
This  normalization  will  make  comparisons  between  signals  and  percent  overlap  more  meaningful. 

Comparison  Study 

A  comparison  study  was  performed  to  analyze  the  effects  of  overlap  with  the  amount  of 
reconstruction  error.  The  ten  signals  that  were  utilized  in  the  experiments  of  Chapter  IV  were  used 
in  this  study  (the  plots  of  the  original  data  are  shown  in  Appendix  A).  Eight  values  of  overlap  were 
selected:  1  point,  5%,  10%,  15%,  20%,  30%,  40%,  and  50%.  The  most  meaningful  representation 
of  the  results  is  a  graph  of  the  error  (Li)  versus  the  percent  overlap.  The  following  plots  are  the 
results  of  implementing  the  Generalized  Malvar  wavelet  and  the  ten  different  signals: 

These  figures  demonstrate  a  definite  correlation  between  the  percentage  of  overlap  and  the 
level  of  reconstruction  error.  More  importantly,  the  plots  show  that  there  is  generally  a  point 
where  the  slope  of  the  error  curve  decreases.  The  reconstructed  signals  have  a  sharp  decrease  in 
the  amount  of  error  through  5  to  15%  overlap  and  then  the  slope  begins  to  level  out. 
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L2 r  rrof  L2  Error 


Figure  5.2.  L2  Error  Versus  Percent  Overlap  (a.  “Wash”,  b.  “Dark”) 


Figure  5.3.  L2  Error  Versus  Percent  Overlap  (a.  “S”,  b.  “She”) 
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Figure  5.4.  L2  Error  Versus  Percent  Overlap  (a.  “Ear”,  b.  Sine) 


Figure  5.5.  L2  Error  Versus  Percent  Overlap  (Two  Low  Freq  Sines) 
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Figure  5.6.  L2  Error  Versus  Percent  Overlap  (a.  Convolution  of  Two  Sines,  b.  Damped  Sine) 


Speech  Signals 

Three  different  sentences  from  the  TIMIT  database  where  transformed  using  three  amounts 
of  overlap  and  three  amounts  of  compression.  The  speech  signals  were  overlapped  by  1  point,  10% 
(26  points),  and  50%  (128  points)  and  compressed  by  80%,  90%,  and  95%.  There  was  a  definite 
decrease  in  the  numerical  error  as  the  amount  of  overlap  was  increased.  The  L2  and  RMS  errors  are 
shown  in  table  5.2.  The  RMS  error,  at  80%  compression,  improved  8.3%  from  one  point  overlap  to 
26  point  overlap  and  the  error  improved  16.3%  when  the  128  point  overlap  is  compared  to  the  one 
point.  All  of  the  sentences  showed  the  same  level  of  improvement;  5-10%  at  26  point  and  15-20% 
at  128  point  overlap. 

The  improvement  in  numerical  error  can  also  be  found  in  the  graphical  display  of  the  signals. 
Figure  5.7a  and  b  plots  the  reconstruction  of  an  utterance  with  one  point  and  128  point  overlaps. 
Differences  are  noted  in  the  low  amplitude  sections  of  the  reconstructions.  One  area  of  the  signal  is 
expanded  in  figure  5.8  to  show  the  differences  in  more  detail.  The  128  point  overlap  demonstrates 
a  better  ability  to  reconstruct  the  ends  of  the  signal  and  wherever  the  signal  has  an  amplitude  less 
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Overlap  ( Compression) 

Li 

RMS 

Sentence  1 

1  pt  (80%) 

9.7447xl0-2 

3.2882x10-* 

26  pt  (80%) 

8. 9384xl0-2 

3.0161x10-* 

128  pt  (80%) 

8.1582x10-2 

2.7529x10-* 

Sentence  2 

1  pt  (90%) 

1.5125xl0-1 

5.1367x10-* 

26  pt  (90%) 

1.3476xl0_1 

4.5767x10-* 

128  pt  (90%) 

1.2046xl0-1 

4.0912x10-* 

Sentence  3 

1  pt  (95%) 

2.9503xl0-1 

6.3541x10-* 

26  pt  (95%) 

2.7649x10-* 

5.9547x10-* 

128  pt  (95%) 

2.5293x10-* 

5.4474x10"* 

Table  5.2.  Reconstruction  Errors  Versus  Amount  of  Overlap  and  Percent  Compression 


than  a  few  hundred.  The  higher  overlaps  gained  more  detail  for  the  signal,  but  did  not  improve 
the  peaks  or  overall  shape  of  the  reconstruction. 

The  improvements  found  in  the  numerical  error  and  in  the  plots  were  not  as  evident  in 
the  audio  reconstruction.  The  increased  overlap  did  reduce  some  of  the  background  noise,  but 
the  increase  did  little  to  actually  improve  the  quality  or  intelligibility  of  the  speech  signed.  The 
intelligibility  of  the  sentence  with  one  point  overlap  was  not  very  different  from  that  of  the  signal 
with  128  point  overlap  because  these  low  level  details  are  not  as  important  to  the  human  ear  as 
the  high  level  peaks  and  the  overall  shape  of  the  speech  signal. 

Summary 

The  exact  amount  of  overlap  which  should  be  used  has  not  been  analytically  determined. 
However,  for  speech  signals  it  was  determined  that  a  one  point  overlap  performs  as  well  as  the 
maximum  amount  of  overlap.  The  percent  overlap  used  for  other  types  of  signals  may  be  more 
important  and  will  depend  upon  the  level  of  reconstruction  desired  and  the  limits  on  the  amount 
of  pre/post  processing  to  be  done.  Figure  5.1  demonstrates  that  there  is  a  real  improvement  to  be 
made  by  increasing  the  overlap.  For  Malvar  wavelets,  the  amount  of  overlap  does  not  effect  the 
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Figure  5.7.  Reconstruction  of  the  Sentence,  “A  Boring  Novel  Is  A  Superb  Sleeping  Pill”,  Using 
(a)  1  Point  Overlap  and  (b)  50%  (128  Point)  Overlap 
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samples 


(b) 

Figure  5.8.  Reconstruction  of  the  Same  Sentence  (From  Sample  7500  to  11000)  Using  (a)  1  Point 
Overlap  and  (b)  50%  (128  Point)  Overlap 
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bit  rate.  The  lapped  transform,  for  uncompressed  data,  sends  the  same  “N*  number  of  samples 
regardless  of  the  overlap  amount.  In  other  words,  the  frame  size,  remains  the  same, 

but  the  amount  of  overlap  with  the  adjacent  frames,  [(j,(j+i),  changes.  There  are  many  instances 
where  the  processing  available,  prior  to  transmission  or  after  reception,  is  quit  acceptable  to  allow 
for  greater  overlap.  Increasing  the  amount  of  overlap  will  increase  the  number  of  computations 
required  prior  to  transmission  and  after  reception.  The  area  that  gives  the  best  return  for  the  least 
number  of  computations  is  the  5  to  10%  overlap  region.  Most  signals  have  improved  significantly  by 
the  10%  mark.  When  the  slope  of  the  curve  levels  out,  the  user  is  spending  more  on  computations 
and  getting  less  of  a  return  (in  the  form  of  a  high  quality  reconstructed  signal). 
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VI.  Real-  Valued  Homomorphic  Filtering 


Introduction 

The  primary  use  of  Homomorphic  filtering  in  the  literature  is  to  perform  deconvolution.  In  this 
way,  a  signal  may  be  separated  into  its  components;  whether  they  be  a  transmitted  signal  convolved 
with  the  transmission  path  system  function  or  the  vocal  tract  system  function  convolved  with  the 
vocal  cords  system  function.  The  most  common  form  of  homomorphic  filtering  is  the  complex 
cepstrum.  But,  as  noted  in  Chapter  II,  the  complex  cepstrum  is  computationally  expensive  to 
utilize  because  of  the  phase  discontinuities. 

The  discrete  cosine  transform  (DCT),  discrete  sine  transform  (DST),  and  the  lapped  orthogo¬ 
nal  transform  (with  a  DST-like  basis,  the  Generalized  Sinusoidal  Malvar  Wavelet  (GMW))  (34),  (18) 
are,  by  definition,  real  functions  (28).  Therefore,  if  the  DCT,  DST  or  GMW  is  used  as  part  of  a 
homomorphic  filter,  we  may  reap  the  benefits  of  deconvolution  without  the  problems  of  phase  dis¬ 
continuities.  The  basis  function  used  with  the  Sinusoidal  Malvar  wavelet  in  Suter  and  Oxley’s 
paper  is  a  real- valued  function  (34): 

AM -^*>  [.(»  +  !)£]  (6.1) 

where  k  is  the  basis  number,  N  is  the  window  length,  and  x  is  the  data  position  within  the  given 
interval.  The  homomorphic  filter  would  be  implemented  as  in  figure  6.1;  where  G(f)  represents 
the  forward  Generalized  Malvar  wavelet  (GMW)  and  G(i)  represents  the  inverse  GMW,  x(n)  is 
the  original  signal  and  x’(n)  is  the  reconstructed  signal,  and  L  shall  be  called  the  homomorphic 
coefficients.  Because  the  resulting  coefficients  are  always  real- valued,  there  is  no  phase  information 
to  track.  Is  it  therefore  less  computationally  expensive  to  use  a  real- valued  basis  function  to 
perform  homomorphic  deconvolution?  The  GMW,  by  definition,  uses  more  than  one  window  to 
calculate  its  coefficients.  However,  all  of  the  signal  must  lie  within  a  single  window  when  performing 
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Figure  6.1.  A  Generalized  Malvar  Wavelet  Homomorphic  System 

deconvolution.  Therefore  it  is  not  really  the  lapped  transform  that  is  being  utilized,  but  its  basis 
function  and  a  single  window. 

Development 

When  using  a  real- valued  function,  the  analysis  of  the  deconvolution  can  become  very  complex 
and  there  are  certain  limitations  on  what  signals  can  actually  be  separated. 

For  example,  a  signal  consisting  of  an  original  signal  and  a  simple  echo  (as  in  figure  6.2)  can 
be  represented  as: 

)  =  «(<)  +  <»«(*  -  r)  (6.2) 

where  a  is  the  echo  amplitude  (0  <  a  <  1)  and  r  is  the  time  difference  between  the  original  signal 
and  the  echo  (r  >  0).  Assume  that  the  signal,  s(t),  damps  out  prior  to  the  arrival  of  the  echo,  so 
that  the  two  signals  do  not  overlap.  The  Fourier  transform  of  this  signal  can  be  represented  simply 
as: 

S«ot(tu)  =  S(w)  +  aF[s(t  -  r)]  (6.3) 

Now  the  Fourier  transform’s  shift  property  (F[s(t  —  r)]  =  S(w)e~*wt)  is  utilized  to  yield: 

=  5(w)[l  +  ae~*wt]  (6.4) 
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However,  for  a  homomorphic  transform  (with  a  different  basis  than  the  Fourier  transform)  to  be 
applied,  we  must  first  know  the  time  shift  property  of  the  transform. 


cm)  = 


Thus  for  f(t  +  r)  we  have: 


G(f(t  +  r)) 


(6.5) 


(6.6) 


Now  let  p  —  t  +  r  so  that  dp  =  dt  and  t  =  p  —  r;  also  using  the  trig  identity  that  sin(a  +  b)  = 
sin(a)  cos(6)  —  cos(a)  sin(6)  yields: 

G(/(t  +  r))  =  y/2/M  J  f(p) sin  [*(*  +  5)^7  cos  [»(fe  +  ^)^  dP  (6-7) 

-s/2/M  J  /(p)cos  | *(k  +  i)^:]  sin  |r(*  +  ^)-^r]  dp  (6.8) 

Now  we  will  let  Ge(f(t))  =  \J2/M  f  /(<)  cos[ir(ifc  +  |)^-]d<  so  that: 

G(f(t  +  r))  =  cos  [*  (k  +  0  G(f(t))  -  sin  [t  (k  +  Ge(/(f ))  (6.9) 

The  time  shift  property  for  a  purely  sinusoidal  basis  is  therefore  a  combination  of  the  original  basis 
and  a  cosine  basis. 

Another  difficulty  in  using  a  sinusoidal  basis  is  that  the  convolution  property  of  the  Fourier 
transform  does  not  hold.  The  convolution  property  for  our  transform  is  given  by: 

G[x(t)  *  y(t)]  =  s/m  [G(x(t))  ■  G(y(t))  -  G'(x(t))  ■  G‘(y(t))]  (6.10) 
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This  means  that  the  signal  must  be  able  to  be  directly  expressed  as  two  or  more  terms  like 
our  example  signal  (sl0,(t)  =  s(t)  +  as(<  —  r))  or  the  signal  must  have  either  even  or  odd  symmetry, 
so  that  one  part  of  the  convolution  property  (equation  6.10)  will  be  equal  to  zero.  If  the  signal  has 
even  symmetry  then  the  Ge(x)  ■  Gc(y)  product  will  equal  zero  and,  likewise,  if  the  input  signal  is 
odd  symmetric  then  the  product  of  G(x)  ■  G(y)  will  be  zero. 

The  use  of  a  real- valued  basis  function  has  some  definite  merit  because  there  is  no  need 
to  perform  phase  unwrapping.  However,  it  is  limited  by  the  type  of  signal  that  can  actually  be 
deconvolved. 


Implementation 

The  homomorphic  filtering  system  was  implemented  with  the  modified  Malvar  wavelet  pro¬ 
gram  used  to  perform  the  homomorphic  filtering  (System  1  using  the  natural  logarithm/exponential) 
in  Chapter  IV.  This  program  accomplishes  all  the  steps  of  the  Sinusoidal  Malvar  wavelet  as  out¬ 
lined  in  Chapter  IV  and  as  proposed  by  Suter  and  Oxley  (34).  Input  data  (x[n])was  created  with 
MATLAB  data  files  and  the  DARPA  TIMIT  Phonetic  speech  database. 


Function 

Input 

Output 

Forward  Transform 

*t»] 

Transform  Coeffs  (a) 

Natural  Log 

o’s 

In  o’s 

Inverse  Transform 

lna’s 

Homomorphic  Coeffs  (L) 

Forward  Transform 

L’s 

Transform  Coeffs  (£) 

Exponential 

C’s 

ec 

Inverse  Transform 

ec 

Output  Signal  (x’[n]) 

Table  6.1.  Homomorphic  Filtering  System  with  Malvar  Wavelet 


Deconvolution:  Deconvolution  is  achieved  by  filtering  the  homomorphic  coefficients  (L). 
Does  the  homomorphic  GMW  system  deconvolve  a  process  into  its  basic  signals  like  the  complex 
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cepstrum?  A  signal  such  as  a  damped  sinusoid  with  an  echo  (as  in  figure  6.2)  is  often  used  to 
demonstrate  the  deconvolution  property  of  the  complex  cepstrum. 


0.5 
0.4 
0.3 

amp  o.2 
0.1 
0 

-0.1 

0  250  500  750  1000  1250  1500  1750  2000 

time(msec) 

Figure  6.2.  Damped  Sinusoid  with  an  Echo 
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The  signal  is  first  transformed  by  the  GMW.  These  wavelet  coefficients  are  plotted  in  fig¬ 
ure  6.3.  The  next  step  is  to  take  the  natural  logarithm  of  the  coefficients  (a).  When  using  the 
power  or  complex  cepstrum,  the  problem  of  an  undefined  logarithm  will  only  occur  when  a  coef¬ 
ficient  is  exactly  0.0.  The  same  procedure  as  outlined  in  Chapter  IV  was  used  to  avoid  any  null 
coefficient  when  using  the  power  or  complex  cepstrum. 

The  inverse  lapped  transform  is  taken  on  the  coefficients  following  the  natural  logarithm. 
This  set  of  values  (plotted  in  figure  6.4)  are  the  homomorphic  coefficients.  The  homomorphic 
coefficients  are  log  amplitude  values  plotted  against  time  (known  as  quefrency  for  cepstral  users). 
These  homomorphic  coefficients  can  now  be  filtered  (liftered)  to  deconvolve  the  signal.  The  first 
filter  implemented  is  an  all-pass  filter  to  demonstrate  that  the  process  is  completely  reversible. 


Reconstruction:  The  homomorphic  coefficients  are  transformed  by  the  following  functions: 
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Figure  6.3.  Malvar  Wavelet  Coefficients 
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Figure  6.4.  Homomorphic  Domain 


Figure  6.5.  Perfectly  Reconstructed  Damped  Sinusoid  with  Echo 

1.  Forward  Sinusoidal  Malvar  Wavelet. 

2.  Exponential  (e£). 

3.  Inverse  Sinusoidal  Malvar  Wavelet. 

The  resulting  signal  is  a  perfect  reconstruction  (with  a  normalized  RMS  error  =  6.3587xlO-10)  of 
the  original  damped  sinusoid  and  echo  (see  figure  6.5). 

In  order  to  remove  the  echo,  the  high  quefrency  information  must  be  removed.  The  desired 
signal  (original  damped  sine)  is  a  high  frequency  signal  that  is  transformed  to  the  low  quefrency 
during  the  homomorphic  filtering  process.  A  low-pass  filter  is  used  on  the  quefrency  domain  to 
eliminate  the  high  quefrency  (low  frequency)  components  of  the  signal.  The  high  frequency  (low 
quefrency)  information  is  maintained  and  the  original  50  Hz  signal  is  reproduced  as  shown  in 
figure  6.6. 
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Figure  6.6.  Reconstructed  Damped  Sinusoid  with  No  Echo 


Summary 

This  chapter  demonstrated  that  it  is  possible  to  utilize  standard  deconvolution  procedures 
with  a  real- valued  homomorphic  system.  The  primary  difficulty  with  this  procedure  is  that  only 
certain  types  of  signals  are  candidates  for  deconvolution.  Signals  that  can  be  deconvolved  with  a 
real-valued  system  are  limited  to  even/odd  symmetric  signals  or  a  signal  that  is  able  to  be  directly 
expressed  as  two  or  more  terms  (««,<(<)  =  s(t)  -)-  os(<  —  r),  where  s(t)  is  sufficiently  damped). 
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VII.  Conclusions 


Summary  and  Conclusions 

The  coding  of  the  algorithms  and  the  experiments  run  with  these  programs  encompass  the 
most  important  contributions  of  this  thesis  effort.  During  the  coding  effort  of  the  Complex  Malvar 
wavelet  algorithm,  the  theory  behind  the  transform  was  validated.  This  untested  algorithm  was 
coded  and  bugs  were  worked  out  of  both  the  program  and  the  theory.  The  validation  effort  was  used 
to  ensure  that  the  program  (listed  in  Appendix  B)  would  lead  to  perfect  reconstruction  of  a  signal. 
Perfect  reconstruction  and  orthonormality  was  accomplished  and  compression  experiments  were 
run  with  the  CMW.  With  the  exception  of  rectangular  windowed  short  time  Fourier  transforms 
(STFT),  this  is  the  first  completely  orthogonal  STFT  to  be  implemented  and  tested  for  digital  signal 
processing  applications.  The  CMW  provides  perfect  reconstruction  of  signals  and  can  reconstruct 
signals  fairly  well  at  high  compression  ratios,  it  was  found  that  the  CMW  can  reconstruct  intelligible 
speech  signals  at  data  rates  below  2400  bps.  These  rates  were  reached  by  compression  of  the 
transform  coefficients,  reducing  the  number  of  bits/coefficient,  and  keeping  positional  information 
to  a  minimum.  The  CMW  bay  be  even  more  appropriate  for  image  and  RADAR  processing. 
Unfortunately,  because  the  CMW  was  defined  so  late  in  the  thesis  cycle,  the  transform  was  not 
extended  to  two  dimensions  or  tested  with  image  processing  applications. 

The  Real- Valued  Malvar  wavelet  software  was  modified  and  improved  to  provide  more  accu¬ 
rate  operation  and  to  incorporate  more  features.  The  original  code  was  not  thoroughly  validated 
by  Raduenz  and  had  a  few  problems  that  did  not  allow  for  perfect  reconstruction.  The  program 
and  theory  were  tested  and  scrutinized  to  identify  and  eliminate  the  bugs.  The  testing  ensured 
that  the  Ada  software  would  provide  an  identical  reconstruction  of  any  type  of  signal  with  varying 
size  windows  and  overlaps.  Previous  research  in  the  area  of  lapped  transforms  has  demonstrated 
perfect  reconstruction  for  data  using  a  50%  overlap,  and  mathematically  derived  conditions  for  a 
more  general  size  overlap.  This  thesis  implemented  a  lapped  transform  algorithm,  based  on  the 
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work  of  Suter  and  Oxley,  that  achieved  perfect  reconstruction  of  a  signal  with  a  one  point  over¬ 
lap  between  intervals.  This  result  is  significant  because  a  transform  based  on  a  one  point  overlap 
requires  half  tie  processing  of  a  transform  based  on  50%  overlap. 

Three  different  improvements  were  made  once  the  basic  program  was  providing  error  free 
operations.  The  basic  program  performed  finite  length  window  processing  by  copying  the  first 
window  of  data  to  the  end  of  the  file  and  then  transforming  an  extra  window  of  data.  The  new 
program  is  able  to  obtain  perfect  reconstruction  of  the  first  and  last  interval  without  copying  data 
or  performing  multiple  calculations  on  one  window.  The  program  “wraps”  the  data  around  so 
that  the  first  and  last  data  point  become  neighbors  and  overlapping  of  the  end  windows  is  done 
directly  between  the  two  windows.  This  is  the  first  technique  that  allows  for  a  finite  length  signal 
to  be  perfectly  reconstructed  using  the  lapped  transform  without  artificially  expanding  the  signal. 
The  Ada  software  for  the  non-expansive  lapped  transform  is  listed  in  Appendix  C.  A  package  was 
written  to  allow  the  Malvar  wavelet  programs  to  recursively  find  the  level  of  thresholding  necessary 
to  provide  a  given  percent  compression.  This  enables  the  researcher  to  qualitatively  compare  the 
results  of  different  systems  at  a  given  level  of  compression.  The  second  Ada  package  written  for 
the  wavelet  programs  provides  the  error  criteria  for  a  reconstruction.  This  also  gives  the  researcher 
a  better  tool  for  studying  reconstruction  techniques. 

Compression  capabilities  of  different  lapped  transform  systems  were  evaluated  and  compared. 
Raduenz  had  found  that  a  homomorphic  system  (using  the  lapped  transform,  the  natural  logarithm, 
and  a  second  lapped  transform)  had  provided  superior  reconstruction  of  a  speech  signal.  The  three 
systems  (as  shown  in  figures  4.3,  4.4,  and  4.5)  were  validated  to  ensure  they  were  working  properly 
at  zero  percent  compression.  The  systems  were  then  evaluated  for  their  reconstructive  abilities  by 
comparing  numerical  errors,  graphical  displays,  and  audio  reconstructions  of  many  different  types  of 
signals  at  various  amounts  of  compression.  In  all  cases  the  most  basic  system  (figure  4.5)  provided 
the  best  reconstruction  of  the  signals.  This  is  important  because  it  allows  for  much  less  pre/post 
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processing  to  gain  the  best  reproduction  possible.  Data  rates  close  to  3000  bps  were  found  to  provide 
reasonable  quality  and  intelligible  speech  reconstruction  for  the  real- valued  Malvar  wavelet. 

The  lapped  transform  theory  allows  for  any  amount  of  overlap  from  one  point  to  50%  of  the 
size  of  the  window,  but  the  theory  does  not  suggest  what  amount  of  overlap  would  be  the  best. 
This  thesis  reports  the  results  of  multiple  experiments  to  compare  and  contrast  the  amount  of 
overlap  to  the  level  of  reconstruction.  The  results  show  that  the  numerical  and  graphical  error 
tends  to  continue  to  decrease  as  the  amount  of  overlap  increases,  although  this  is  not  the  case  for 
all  signals.  The  rate  at  which  the  error  decreases  typically  leveled  out  between  5  and  15%  overlap. 
This  would  suggest  that  about  10%  overlap  would  be  an  optimal  amount  for  most  signals.  Speech 
signals  showed  this  same  pattern  of  improvement,  however  the  decrease  in  error  could  not  be  heard 
in  the  reconstructed  speech.  The  speech  signal  is  improved  in  low  amplitude  sections  that  are  not 
easily  differentiated  by  the  human  ear.  This  again  suggests  that  the  most  basic  method,  the  lapped 
transform  with  one  point  overlap,  be  used  for  speech  compression  purposes.  This  will  save  half  of 
the  processing  required  for  a  50%  overlap.  Deconvolution  was  also  studied  using  the  basis  from  the 
real-valued  Malvar  wavelet.  Deconvolution  was  shown  to  be  possible  using  a  sinusoidal  function, 
however  there  are  certain  limitations  on  the  signal.  The  signal  must  have  even  or  odd  symmetry  or 
be  able  to  be  directly  expressed  as  two  individual  terms  (such  as  a  damped  signal  with  an  echo). 

Follow-On  Research 

The  extension  of  the  Complex  Malvar  wavelet  to  two  dimensions  is  an  intriguing  possibility 
for  follow-on  work.  The  orthogonality  of  the  CMW  could  provide  excellent  results  with  image 
processing  applications.  The  CMW  could  be  substituted  into  almost  any  algorithm  or  application 
which  utilizes  the  complex  spectrum  of  a  signal.  A  two  dimensional  Ada  program  would  not 
be  trivial  because  the  FFT  (see  Raduenz  (26))  routine  would  also  have  to  be  extended  to  two 
dimensions.  The  CMW  also  has  some  important  one  dimensional  applications  that  should  be 
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researched.  Further  speech  compression  experiments  could  be  performed  using  more  elaborate 
compression  and  coding  techniques.  Data  rates  should  be  calculated,  taking  into  account  the 
quality  of  the  reconstruction,  and  the  results  compared  back  to  a  baseline  system  like  LPC-10.  The 
CMW  can  also  be  utilized  in  RADAR  signal  analysis  and  to  perform  complex  cepstral  processing. 
Deconvolution  or  other  complex  cepstral  applications  may  have  higher  performance  with  the  CMW. 
Besides  the  code  presented  in  this  thesis,  a  phase  unwrapping  package  would  have  to  be  produced 
in  order  to  use  the  CMW  in  a  homomorphic  filter.  Similar  overlap  experiments  could  be  performed 
with  the  CMW.  It  was  found  that  the  overlap  has  little  effect  on  speech  signals,  but  that  may  not 
hold  true  for  images. 

A  hybrid  wavelet  approach  looks  promising  for  future  research.  Wavelet  packets  and  the 
Malvar  wavelet  could  be  used  together  for  digital  signal  processing.  Wavelet  packets  could  be 
used  to  segment  a  data  stream  while  the  Malvar  wavelet  would  then  analyze  the  data  within  each 
segment. 

The  basis  functions,  window  size,  and  amount  of  overlap  are  all  variable  entities  within  the 
Generalized  Malvar  wavelet  theory.  All  work  to  date  has  been  accomplished  using  a  set  value  for 
these  entities  throughout  a  given  signal.  It  may  be  possible  to  vary  the  basis  function,  size  of 
the  window,  and  the  amount  of  overlap  to  the  changes  within  the  signal.  As  the  statistics  of  the 
signal  are  changing,  different  basis  functions  may  yield  better  reconstruction  than  others.  Further 
compression  may  also  be  obtained  by  scaling  the  size  of  the  window  to  the  type  of  signal.  A  signal 
that  is  not  changing  very  rapidly  could  utilize  a  large  window,  while  a  signal  that  changes  very 
quickly  would  probably  be  best  transformed  with  a  smaller  window. 

A  comparison  of  the  lapped  transform  with  a  speech  coder  like  LPC-10  would  provide  greater 
insight  into  the  capabilities  of  the  algorithm.  There  are  many  advanced  coding  techniques  that 
were  not  utilized  with  the  lapped  transform  in  this  thesis.  If  these  advanced  coding  techniques 
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could  be  incorporated  with  the  wavelet  compression,  then  a  better  comparison  could  be  made.  The 
wavelet  technique  may  be  able  to  provide  superior  quality  and  intelligibility  for  less  than  2400  bps. 
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Figure  A.2.  Utterance  of  the  Word  (a)  "  Vash”  and  (b)  “She” 


Figure  A.4.  (a) 


Figure  A.6.  (a)  Convolution  of  Two  Different  Frequency  Sinusoids  and  (b)  Damped  Sinusoid  with 

an  Echo 
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Appendix  B.  Ada  Source  Code  For  The  Complex  Valued  Afalvar  Wavelet 


This  appendix  provides  the  source  code  used  during  the  development  and  testing  of  the 
Complex  Valued  Malvar  wavelet. 


- 10 - 20 

with  Text.Io; 
with  Complex.Pkg; 
use  Complex.Pkg ; 
with  Math.Lib; 
use  Math.Lib; 
with  Vector .Package; 
use  Vector .Package ; 
with  Type.Package ; 
use  Type.Package ; 
with  Print .Package; 
use  Print .Package; 
with  FFT.Pack; 
use  FFT.Pack; 


procedure  caw  is 

—  Performs  a  forward  complex  Malvar  Vavelet 

—  and  an  Inverse  complex  Malvar  wavelet. 

—  Inputs  must  be  real,  transform  coefficients  are  real  and  imaginary. 

—  Output  values  are  real. 

package  Integer.Io  is  new  Text.Io. Integer.Io(integer); 
package  Float.Io  is  new  Text.Io. Float.Ioff loat) ; 

Max.Data.Length  :  constant  :=  100.000; 

Max_Partitions  :  constant  :=  500; 

Big_Daddy_Data_Vector  :  Real.Vector  (1 . .Max.Data.Length) ; 

Big_Daddy_Partition_Vector  :  Partition.Vector  (0. .Max.Partitions) ; 

Actual.Data.Length  :  integer  :=  0; 

Actual.Partitions  :  integer  :=  0; 

Vindow.Size  :  integer  :=  0; 

Over lap. Amount  :  integer  :=  0; 

Alpha.Iame , 

Outfile.Iame  :  string  (1..30); 

Alpha.lame.Length , 

Outf ile.Iame.Length  :  integer  range  1..30; 


- 30 - 40 - 50 - 60 - -70 - 78 

78* 
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Get.File.Data  prompts  the  user  lor  the  prepared  info  file  name 
and  then  reads  in  the  data  in  the  foliosing  order 

(1)  Input  Data  Filename 

(2)  Derivative  Output  Filename 

(3)  Alpha  Output  Filename 

(4)  Output  Data  Filename 

(6)  Input  Output  Difference  Filename 

(6)  General  Info  Filename 

(7)  lumber  of  Points  in  the  Data  Vector 

(8)  lumber  of  required  partitions 

(9)  -  (?)  loop  for  input  given  in  (8)  and  get 
Boundary  point  then  Overlap  for  each  partition 


procedure  Get_File_Data 


(Alpha.Filename 
Alpha_Length 
Output .Filename 
Outname.Length 
Partitions 
Windos.Size 
Qverlap.Amount 
lumber.Of .Part it ions 
Data 
Points 


in  out  string; 

in  out  integer; 

in  out  string; 

in  out  integer; 

in  out  Partition.Vector; 

in  out  integer; 

in  out  integer; 

in  out  integer; 

in  out  Real .Vector; 

in  out  integer  )  is 


Filenames,  lengths  must  be  in  out  so  they  can 
be  written  to  info  file. 

lumber.Of.Partitiona  and  Points  must  be  in  out 
Parameters  because  they  are  used  in  a  loop  to 
Read  in  the  appropriate  amount  of  data 

Infile, 

Datafile 
Temp. Input 
Data.Filename 
DataJFilename.Length 
In.Filename 
InJFilename.Length 
Counter 

begin 

Text.Io . lew.Line ; 

Text.Io. put ("What  is  the  name  of  the  file  containing  "); 
Text.Io. put ("the  prepared  information  ?  >  "); 

Text.Io . Get .Line (In.Filename ,  In.Filename.Length) ; 

Text.Io . Open (Inf ile ,  Text.Io . In.File, 

In.Filename ( 1 . . In.Filename.Length) ) ; 


lext.io .  r  ne_  type ; 
integer; 
string(l. .30) ; 
integer; 
string(l. .30) ; 
integer; 
inteser  :=  0: 
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—  low  I  call  in  the  required  user  data 

Text.Io. Get_Line(Inlile,  Data.Filename,  Data.Filename.Length) ; 
Text.Io. Get_Line(Inlile,  Alpha.Filename,  Alpha.Length) ; 
Text.Io. Get_Line(Inlile,  Output .Filename,  Out name .Length) ; 
Integer.Io. get (Inf ile, Points) ; 

Integer.Io . get ( Ini ile , Window.Size) ; 

Integer _Io . get ( Ini ile , Over lap. Amount ) ; 


Iumber.01 .Part it ions  :=  Points/Window.Size; 

lor  j  in  1 . .lumber.Ol.Partitions  loop 
Counter  :=  Counter  +  1; 

Partitions(j) .Boundary  :=  Window_Size*Counter; 
Partitions(j) .Overlap  :=  Over lap. Amount; 
end  loop; 

Text.Io . Close (Ini ile) ; 

Part it ions (0) .Boundary  :  =  1; 

Partitions (0) .Overlap  Overlap.Amount ; 


—  now  I  call  in  the  numbers  lor  the  vector  Data 
Text.Io . Open(Datal ile .Text.Io . In.File , 

Data.Filename (1 . . Data. Filename. Length) ) ; 

lor  j  in  1.. Points  loop 

Float_Io.Get(Datalile,  Data(j));  —  lor  lloating  point  inputs 
— Integer.Io .Get(Datalile .Temp.Input) ; 

— Data(j)  :=  lloat (Temp.Input) ;  —  lor  integer  inputs 

end  loop; 

Text.Io.Close (Data! ile) ; 


end  Get.File.Data; 

- 10 - 20 - 30 - 40 - 50 
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—  Mult iply.By .Window 


(Data 

:  in  Real. Vector; 

Folded.Data 

:  in  out  Real. Vector; 

Window. j 

:  in  out  Real.Vector; 

A-j 

:  in  integer; 

A.j.plus.l 

:  in  integer; 

Epsilon 

:  in  out  integer; 

Window.size 

:  in  integer; 

J 

:  in  integer)  is 

Pi 


:  constant  :=  3.141592654; 
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begin 


Window.Boundary 
Start .File 
End.File 
Counter 
n_plua_l 


constant  float  :=  1.0/sqrt(2.0) ; 

:  integer  :=  Data '1  .'•at; 

:  integer  .=  Data 'last; 

:  integer  :=  0; 

:  integer  :  =  window. size ; 


- *****Hultiply  By  Window***** - 

—  Left  edge 

if  j  *  0  then  — A.  j  =  0 

Folded.DataU.  j )  :=  2.0  *  Data(End.File)  *  Windon_j(0); 
counter  :=  1; 

for  x  in  (A_j+1) . • (A.j  +  EPSIL0I)  loop 
Folded.Data(x)  :=  Data(x)  *  Window. j (counter)  + 
Data(End_f ile  -  x)  *  Window.j (-counter) ; 
counter  :  =  counter  +  1; 
end  loop; 
else 

counter  : =  0; 

for  x  in  (A_j)..(A_j  +  EPSIL0I)  loop 
Folded.Data(x)  :=  Data(x)  *  Window.j (counter)  + 

Data(2*A_j  -  x)  *  Window. j (-counter) ; 
counter  :*  counter  +  1; 
end  loop; 
end  if; 


—  Center 

for  x  in  (A_j  +  EPSILOI  +  1) . . (A.j.plus.l  -  EPSILOI  -  1)  loop 
Folded.Data(x)  :=  Data(x); 
end  loop; 

—  Right  edge 

if  A.j.plus.l  =  Data’last  then  --(j  =  partitions ’last-1) 
counter  :  =  EPSILOI; 

for  x  in  (A.j.plus.l  -  EPSILOI) .. (A.j.plus.l  -  1)  loop 

Folded.Data(x)  :=  Data(x)  *  Window.j (n.plus.l  -  counter)  - 
Data(0  +  counter)  * 

Window.j (n_plus_l  +  counter); 
counter  :=  counter  -  1; 
end  loop; 

Folded_Data(A_ j .plus.l)  :=  0.0; 
else 

counter  :=  EPSILOI; 

for  x  in  (A.j.plus.l  -  EPSILOI) .. (A.j.plus.l)  loop 

Folded.Data(x)  :=  Data(x)  *  Window.j (n.plus.l  -  counter)  - 

Data(2*A_j_plus_l  -  x)  *  Window.j (n.plus.l  +  counter); 
counter  :=  counter  -  1; 
end  loop; 
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end  if; 


end  Multiply.By.Window; 

- 10 - 20 - 30 - 40 - BO - 60 - 70 - 78 

—  Compute.Coeff icients 

—  This  procedure  takes  in  data  from  0  -  1-1  or  A.j  -  A_j+1-1  (I  points)  and 

—  returns  coefficients  stored  from  0  -  1-1  or  A.j  -  (A_j+1  -  1)  (■  points) 


procedure  Compute.Coefficients  (Data_Segment  :  in  Real.Vector; 

FFT_0ut  :  out  Complex. Vector; 

I  :  in  integer)  is 

Beta  :  Real.Vector (Data.segment1 range) ; 

Temp.FFT  :  Complex.Vector (Data.segment ’range) ; 

—  Data.Segment1 range  =  (Start .Data) .. (End.Data  -  1) 

begin 

—  (1)  Scale  Data_Segment  Prior  to  taking  FFT 

Beta(Beta'f irst)  : =  Data.segment (Data.segment ’first) ; 
for  x  in  l..(I-l)  loop 

Beta(Beta1fir8t  +  x)  := 

2 . 0*Data.Segment (Data.segment 1 first  +  x); 

end  loop; 

—  (2)  Make  Data.segment  Complex  Valued 
Temp.FFT  : =  Complex.Of ( (Beta) ) ; 

—  (3)  Perform  Inverse  FFT  (See  FFT.Pack) 

FFT  (Temp.FFT,  True) ;  —  includes  1/1  scaling 


—  (4)  Scale  data  by  1/2  (Inv  FFT  includes  1/B  scaling,  but 
we  require  1/21) . 

for  x  in  FFT.Out 'range  loop  — (Start.Data) . . (End_Data-l) 
FFT.Out(x)  :=  0.5*Temp_FFT(x) ; 
end  loop; 

end  Compute.Coeff icients; 
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procedure  Reconstruct ion_FFT(FFT_In  in  Complex.Vector; 
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FFT_Out  out  Co*plex_Vector)  is 

Tenp_FFT  :  Coaplex_Vector(FFT_In* range)  :  =  FFT_In; 
begin 

—  (1)  Parlors  FFT  (See  FFT.Pack) 

FFT(Tenp_FFT,  False); 

lor  x  in  FFT.Out 'range  loop 
FFT_Out (x)  :=  Tenp_FFT(x) ; 
end  loop; 

—  Scale  Iron  Beta  values  to  H_nl  values. 

FFT_out(FFT_Out* first)  :=  2 . 0*Temp_FFT(Tenp_FFT’lirst) ; 

end  Reconstruction_FFT; 


—  Divide_By_Windovs 


procedure  Divide_By_ Windows  (  TenpJJata 

Data_Segmant 

Window 

A_j_plus_l 

EPSILOI 

Window.Size 

J 


:  in  Complex.Vector; 

:  out  Conplex_Vector; 
:  in  Real.Vector; 

:  in  integer; 

:  in  integer; 

:  in  integer; 

:  in  integer; 

:  in  integer)  is 


Counter 

End_File 

n_plus_l 

test 


:  integer  :=  0; 

:  integer  :=  Temp_Data'last; 
:  integer  :=  window.size; 

:  float; 


begin 


—  Calculate  Data_Segaent(A_j  -  e_j  +1..  A_j  -  1) 
if  j  =  0  then  — A_j  *  0 

counter  : =  EPSILOI; 

for  i  in  (A_j  -  EPSILOI  +  l)..(A_j  -  i)  loop 

Data_Segment(i)  :=  Window (n_plus_i  -  counter  +  1)  * 
Temp_Data(End_File  +  i) 

+  Window(n_plus_l  +  counter  -  1)  * 
Temp_Data(2*A_j  -  i) ; 
counter  :=  counter  -  1; 
end  loop; 
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else 


counter  :=  EPSIL01; 

for  i  in  (A.j  -  EPSILOI  +  l)..(A_j  -  1)  loop 

Data_Segment(i)  : =  Window(n_plus_l  -  counter  +  1)  * 
Tenp_Data(i)  + 

Window (n.plus.l  +  counter  -  1)  * 
Temp_Data(2*A_j  -  i); 
counter  :=  counter  -  1; 
end  loop; 
end  il; 


—  Calculate  Data_Segment(A_j) 
test  :=  2.0  *  window(O); 

Data_Segnent (A. j )  :=  Teap_Data(A_j)/test; 

—  Temp.Data  was  calculated  Iron  A_(j)-e  to  A_(j+1)-1,  so  there 
—  is  no  Temp_Data(End_File) 

—  Calculate  Data_Segaent (A_ j  +  1  ..  A_j  +  e_j  -  1) 
if  j  =  0  then  — A_j  =  0 

counter  : =  1; 

for  i  in  (A_j  +  l)..(A_j  +  EPSILOI  -  1)  loop 

Data.Segment(i)  :=  tfindow(counter)  *  Teap_Data(i)  - 

Hindow(-counter)  *  Temp_Data(End_File  -  i); 
counter  :=  counter  ♦  3; 
end  loop; 
else 

counter  :=  1; 

for  i  in  (A.j  +  1) . . (A_j  +  EPSILOI  -  1)  loop 

Data.Segment(i)  :=  Window (counter)  *  Temp_Data(i)  - 
Window (-counter)  *  Temp.Data (2*A_j  -  i); 
counter  :=  counter  +  1; 
end  loop; 
end  if; 


—  Calculate  Data_Segment (A_ j  +  e_j  ..  A_j+1  -  e_j) 
for  i  in  (A_j  +  EPSILOI) . . (A_j_plus_l  -  EPSILOI)  loop 
Data.Segnent(i)  :=  Tenp.Data(i) ; 
end  loop; 

end  Divide.By.Windows; 


procedure  Error.Out  (Recon.Oata  :  in  Complex. Vector; 

Actual.Data  :  in  Real.Vector; 

L2.Error  :  in  out  float; 

RNS.Error  :  in  out  float; 

L.IIF.Error  :  in  out  float; 
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lora.RIV  :  in  out  float)  i> 


error. 

sual. 

sua2. 

sob3. 

count 

:  float  :*  0.0; 

Recon.RIV 

:  float  :=  0.0; 

Actual.RIV 

:  float  :=  0.0; 

begin 

for  1  in  Actual.Data* range  loop 

if  abs(Actnal_Data(l)  -  ReconJ)ata(l) .real)  >  error  then 
error  :=  abs(Actnal_Data(l)  -  Recon_Data(l) .real) ; 
end  if; 

stuil  :=  Actual_Data(l)*Actual_Data(l)  +  stuel; 

suu2  :  =  Recon_Data(l)  ,real*Recon_Data(l)  .real  +  sun2; 

sum3  :=  (Actnal_Data(l)-Recon_Data(l) .real)* 

(Actnal_Data(l)-Recon_Data(l) .real)  +  sun3; 
count  :=  count  +  1.0; 
if  1  =  1  then 
suml  :=  0.0; 
sum2  :=  0.0; 
sui>3  :=  0.0; 
error:*  0.0; 
count  :*  0.0; 
end  if; 

end  loop; 

RMS_Error  :*  sqrt(sum3) /count; 

L2_Error  :=  sqrt(&bs(suiBl-sun2)/sual) ; 

L_IIF_Error  :=  error; 

Recon_RIV  :=  ln(sua2)/ln(2.0) ; 

Actual_RIV  : =  ln(sual)/ln(2.0) ; 

■orm.RIV  :=  Recon_RIV/Acttxal_RIV; 

end  Error _0ut; 
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—  Do_Work 

procedure  Do.Vork  (  Data 

in  Real .Vector; 

Partitions 

in  Partition.Vector; 

Vindov.Size 

in  out  integer; 

Epsilon 

in  out  integer; 

Alphaf ile.String 

in  String; 

Outf ile.String 

in  String  )  is 
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Pi 

Transforaed.Data 

Tenp.data 

Window 

FFT.Data, 

Output .Data 
Multiply 


:  constant  :=  3.141692664; 

Real.VectorC (Data’f irst-1) . . (Data 'last)) ; 

Conplez.Vector 

( (Data’f irst-1) . .Data'last); 

Real.Vsctor 

((-EPSILOI) . . (Vindov.Siza  +  EPSILOI)) ; 

Conplsz.Vactor  ((Data’f irst-Partitions(O) .Owerlap-1) . . (Data'l 
:  Raal_Vsctor( (Data’f irst-1). . (Data'Last-1)) ; 


Start .Data,  End.Data, 
Start .Window,  End. Window, 
Last.Output .Point , 
Last.Window.Point , 
Counter,  J.Count ,  Temp 


—  A.j  and  A_j+1 

—  A.j  -  a.j  and  A_j+1  +  e_j+l 

—  A.j+1  -  a_j+l 

*-j  +  •-) 

:  integer  :=  0; 


Alphafile, 

Outfile 

L2_Error 

L.lIF.Error 

RMS.Error 

lorm.RIV 

Percent.Compression 

Threshold 

Percent.Conpressed 

CompresBion.Count 

lo.Zeros 

■o.Zeros.R 

Maz 

Adjust 

TooHigh 

TooLow 

Start .File 

End.File 

n 

n_plus_l 


Tezt.Io . File.Type ; 

float  :=  0.0; 

float  :  =  0.0; 

float  :=  0.0; 

float  :=  0.0; 

float  :  =  0.0; 

float  :=  0.0; 

float  :  =  0.0; 

integer  :=  0; 

integer  :  =  0; 

float  : =  0.0; 

float  : =  0.0; 

float  : =  1.0; 

integer  :=  0; 

integer  :=  0; 

integer  : =  Data 'first; 

integer  :  =  Data'last; 

constant  integer  :  =  0; 

constant  integer  : =  Window.Size; 


begin 

Tezt.Io . Create  (Alphafile,  Tezt.Io. Out .File,  Alphaf ile.String) ; 
f«xt_Io . Create  (Outfile,  Tezt.Io. Out .File,  Outf ile.String  ); 


Tezt.Io . lew.Line ; 

Tezt.Io. put ("What  is  the  desired  percent  conpression  ?  "); 
Tezt.Io . leu. line ; 

Tezt.Io. put ("(ie  0.86  for  85%  conpression)”); 

Tezt.Io . lew.line ; 

Float _Io . get (Percent .Conpression) ; 

Tezt.Io . lew.line ; 
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♦♦♦♦♦CREATE  WIIDOW***** - 


—  left  rising  sdgs 
if  EPSILOI  ■  0  then 

Vindos(n)  :=  1.0/(sqrt(2.0)); 

•Is* 

for  x  in  (n  -  EPSILOI).. (n  +  EPSILOI  -  1)  loop 
Window(x)  : =  sin((Pi/(4.0*float(EPSIL0I)))  ♦ 

(float (x  -  n  +  EPSILOI))); 

snd  loop; 
and  if; 

—  cantor 

for  x  in  (n  +  EPSILOI) . . (a.plus.l  -  EPSILOI)  loop 
Vindow(x)  :=  1.0; 
and  loop; 

—  right  falling  adga 
if  EPSILOI  =  0  than 

Window(n.plua.l)  : =  1.0/(sqrt(2.0)); 

alsa 

for  x  in  (n_plns_l  -  EPSILOI  +  1) . . (a_plus_l  +  EPSILOI  - 
loop 

Vindow(x)  :=  cos((Pi/(4.0efloat(EPSIL0I)))  ♦ 

(float (x  -  n.plus.l  +  EPSILOI))); 

and  loop; 

Vindow(n_plus_l  ♦  EPSILOI)  :=  0.0; 
and  if; 


J_Connt  :=  -1; 

for  j  in  0. . (Partitions'laat-l)  loop 
J_Connt  :=  J_Connt  +  1; 

begin  —  each  segment  is  transformed  within  this  loop 
if  j .count  =  0  then 
TEMP  :=  1; 

alsa 

TEMP  :=  0; 
and  if; 

- Assign  j  Dependant  Values - 

Start.Data  :  =  Part  it  ions(j)  .Boundary  -  Taap; 

End_Data  :=  Partitions( j+1) . Boundary; 

Start .Window  : =  Start .Data  -  EPSILOI; 

End .Window  :=  End .data  +  EPSILOI; 

- 10 - 20 - 30 - 40 - 50 - 60 - TO - 

Text _Io. put ("Working. . . .  Processing  Data  Segment"); 
Integer_Io.put( j+1) ; 

Taxt.Io . lew.Line ; 


—78 
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Mult iply_By_Window(Data(Start .File. .End.File) , 

Transforaed_Data( (Start .Data) . .End.Data) , 
Window (window’ range) ,  Start  .data, 
End.Data,  EPSILOI,  Window.Siza,  J.Count); 


Coaputa.Coaff icients (Trans! ora*d_Data( (Start.Data) . . (End_Data-l)) , 

FFT_Data((Start_Data) . . (End.Data-1)) , 
Window.Siza) ; 

—  This  rout in#  performs 

—  (1)  acala  th«  data, 

—  (2)  an  invars*  FFT, 

—  (3)  assigns  coefficients  0  -  1-1  to  Transforaed_Data 

- Print  Alpha  Filas - 

lor  x  in  (Data’lirst-1) . . (Data’ last-1)  loop 

—  Print  output  valuas  to  tha  alpha  lila  (coalliciants) . 
Float.Io.Put  (Alphafila,  FFT.Data(x) .raal) ; 

Text.Io. Ies_Line( Alpha! ila) ; 
end  loop; 

10 - 20 - 30 - 40 - 60 - 60 - 70 - 78 

-  REC0ISTRUCTI0I  - 


Reconstruct ion.FFT  (FFT_Data(Start_Data. . (End_Data-l) ) , 

Taap_Data(Start_Data. . (End_Data-l))) ; 
—  This  routine  parforas  a  forward  FFT 

end; - block 

and  loop;  —  —  j  loop 

J.Count  : =  -1 ; 

lor  j  in  0. . (Partitions ’last-1)  loop 
J.Count  :=  J.Count  +  1; 

begin  —  each  sagaant  is  transforaed  within  this  loop 
if  j .count  =  0  than 
TEMP  :=  1; 

else 

TEMP  :=  0; 
and  il; 

- Assign  j  Dependant  Valuas - 

Start.Data  : =  Partitions(j) .Boundary  -  Taap; 

End.Data  : =  Partitions(j+l) .Boundary; 

Start .Window  Start.Data  -  EPSILOI; 

End .Window  : =  End.Data  +  EPSILOI; 
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Laat .Output .Point  :=  End.Data  -  EPSILOI; 

Last.Vindou.Point  : =  Start .Data  +  EPSILOI; 

- 10 - 20 - 30 - -40 - 50 - 60 - TO - 

Diride.By. Windows (Tesp_Data(Taaq>.Data 'rang*) , 

Output .Data (Start .Window. .Laat. Output .Point) , 
Window (window ’range) , 

Start .Data,  End.Data,  EPSILOI, 

Vindow.Siza,  J.Count); 

if  j. count  =  0  than 

for  k  in  0.. (EPSILOI  -  1)  loop 

Output_Data(End_fila  -  k)  :=  Output .Data(-k) ; 
and  loop; 
and  if; 


and; - block 

and  loop;  - j  loop 


- CALCULATE  OUTPUT  ERRORS  - 

Error.Out  (  Output .Data(Data’range), 

Data (Data’ range) , 

L2_Error,  RMS.Error,  L.IIF.E rror,  lom.RIV); 


- PRUT  OUT  TO  FILES - 

Taxt.Io . Close (Alphaf ila) ; 

Taxt.io. Put .Lina ("print  output  to  file"); 

Taxt.Io . law.Lina ; 

Taxt.io. put ("L2  Error  =  "); 
Float_Io.put(L2_Error) ; 

Taxt.io . law.Lina ; 

Taxt.io. put("L(Inf inity)  Error  *  "); 

Float.Io.put(L.IIF.Error) ; 

Taxt.io . law.lina ; 

Taxt.io. put ("RMS  Error  =  "); 
Float.Io.put(RMS.Error) ; 

Taxt.io . law.Line ; 

Taxt.io. put ("lornalized  RIT  =  "); 

Float.Io . put (lora.RIV) ; 

Taxt.io . law.Line ; 

for  x  in  Data 'range  loop 

Float.Io. put(0utf ila, Output.Data(x) .real) ; 

—  Real  signals  in  — >  Real  signal  out. 

- Float.Io .put (Outf ila .Output .Data(x) .  iaag) ; 

Taxt.io. Iev_Line(Outf ila) ; 
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•ad  loop; 


Text_Io.Close(Oatf ile) ; 
•ad  Do.Vork; 


begin  — aaia 

Get.File.Data  (Alpha.Iaae ,  Alpha. lame.Length, 

Oatlil«_Iaa«,  Outf ile.Iaae.Length, 
Big_Daddy_Partition_Vector,  Vindov.Size, 
Overlap.Amount ,  Actual.Partitions, 
Big_Daddy_Data_Vector ,  Actual.Data.Length) ; 

—  By  sanding  th«  correct  size  arrays  into  Do.Work,  tha  array 

—  attributes  caa  be  used  to  determine  the  size  rather  than 

—  passing  another  parameter 

Do. Work  (Big_Daddy_Oata_ Vector ( 1 • . Actual.Data.Length) , 

Big_Daddy_Partition_Vactor  (0. .ActualJPartitions) , 
Vindov.Size,  Overlap.Amount , 

Alpha. lame  ( l . .  Alpha_lame.Length) , 

Outf ile_Ia*e(l . .Outf ile.Iaae.Length)) ; 

end;  — main 
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Appendix  C.  Ada  Source  Code  For  The  Real  Valued  Malvar  Wavelet 


This  appendix  provides  the  source  code  used  during  the  development  and  testing  of  the  Real 
Valued  Malvar  wavelet. 


- 10 - 20 - 30 - -40 - 50 - 60 - 70 - 78 

— Iotas  :  - 

with  Taxt_Io;  —  _____  73* 

with  Complex.Pkg; 

use  Complex.Pkg; 

with  Hath_Lib; 

use  Hath_Lib; 

with  Vector .Package ; 

use  Vector.Package; 

with  Type.Package ; 

use  Type.Package; 

with  Print .Package ; 

use  Print_Package ; 

with  FFT_Pack ; 

use  FFT_Pack; 


procedure  glot  is 

—  Performs  the  basic  GLOT  and  Inverse  GLOT 

—  Performs  the  overlapping  without  copying  the  end  point  data. 

—  The  program  "wraps"  the  data  around,  so  that  the  first  data 

—  point  connects  with  the  last  data  point. 

package  Integer.Io  is  new  Text_Io.Integer_Io(integer) ; 
package  Float.Io  is  new  Text_Io.Float_Io(float) ; 

Max.Data.Length  :  constant  :  =  100.000; 

Nax.Partitions  :  constant  :=  500; 

Big_Daddy .Data. Vector  :  Real.Vector  (1. . Max_Data_Length) ; 

Big.Daddy .Partition. Vector  :  Partition.Vector(0. .Nax.Partitions) ; 

Actual.Data.Length  :  integer  :=  0; 

Actual.Partitions  :  integer  :  =  0; 

Alpha.Iame , 

Outfile.Iame  :  string  (1..30); 

Alpha.Iame .Length , 

Outfile.Iame .Length  :  integer  range  1..30; 
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Get.File.Data  prompt*  the  user  for  the  prepared  info  file  name 
and  then  reads  in  the  data  in  the  following  order 

(1)  Input  Data  Filename 

(2)  Alpha  Output  Filename 

(3)  Output  Data  Filename 

(4)  lumber  of  Points  in  the  Data  Vector 

(5)  lumber  of  required  partitions 

(6)  -  (?)  loop  for  input  given  in  (S)  and  get 
Boundary  point  then  Overlap  for  each  partition 


(Alpha.Filename 

:  in  out 

string; 

Alpha.Length 

:  in  out 

integer; 

Output .Filename 

:  in  out 

string ; 

Outname.Length 

:  in  out 

integer; 

Partitions 

:  in  out 

Partition.Vector; 

lumber _0f .Part it ions 

:  in  out 

integer; 

Data 

:  in  out 

Real.Vector; 

Points 

:  in  out 

integer  )  is 

Filen=”ne  ,  lengths  must  be  "in  out"  so  they  can 
be  written  to  info  file  (if  want  to  use  an 
inf irmation  file. 

lumber _0f .Partitions  and  Points  must  be  "in  out" 
parameters  because  they  are  used  in  a  loop  to 
read  in  the  appropriate  amount  of  data. 

Infile, 

Datafile 
Temp.Input 
Data.Filename 
Data_Filename_Length 
In.Filename 
In_Filename_Length 
Counter 

begin 

Text.Io . lew.Line ; 

Text.Io. put ("What  is  the  name  of  the  file  containing  "); 
Text_Io.put("the  prepared  information  ?  >  "); 
Text_Io.Get_Line(In_Filename,  In_Filename_Length) ; 

Text.Io. Open (Infile,  Text_Io.In_File, 

In.Filename ( 1 . . In_Filename_Length) ) ; 

—  low  I  call  in  the  required  user  data 

Text_Io.Get_Line(Infile,  Data.Filename ,  Data_Filename_Length) ; 
Text_Io.Get_Line(Infile,  Alpha.Filename ,  Alpha.Length) ; 


i ext_io . r lie. xype ; 
integer ; 
stringd .  .30) ; 
integer; 
string(l. .30) ; 
integer; 
integer  :=  0: 
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Text.Io. Get .Line (Inlile,  Output .Filenaae,  Outname.Length) ; 
Integer _Io . get (Ini ile , Points ) ; 

Int  eger _ Io . get ( Ini ile , lumber .01 .Part it ions ) ; 


-  Read  In  Data  Fro*  Data  File  - 

—  Set  up  all  ol  the  Partitions  and  the  amount  ol 

—  overlap  lor  each  window. 

lor  j  in  1. .luaber .Ol.Partitions  loop 

Integer_Io.get(Inlile, Part i t i ons ( j ) . Boundary ) ; 

Integer.Io . get (Inlile , Partitions (j ) . Overlap) ; 
end  loop; 

Text.Io. Close (Inlile) ; 

—  A(0)  position  initially  set  to  1,  later  set  to  0 

—  so  that  there  are  0  ->  I  data  positions. 

Partitions (0) .Boundary  :=  1; 

Partitions(O) .Overlap  : =  Part itions (luaber .ol.Partitions) .Overlap; 

—  Call  in  the  data  lor  the  vector  Data 
Text.Io . Open (Data! ile , Text.Io . In.File , 

Data.Filenaae(l. . Data_Filenaae .Length) ) ; 
lor  j  in  1 . . Points  loop 

Float_Io.Get(Datalile,  Data(j)); 

— Integer.Io . Get (Datalile .Temp.Input) ; 

— Data(j)  : =  lloat (Temp. Input) ;  —  lor  integer  inputs 

end  loop; 

Text.Io. Close(Datal ile) ; 


end  Get_File_Data; 
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—  Multiply .By .Window 


procedure  Multiply .By .Window  (Data 

Folded.Data 

Window.j 

A-j 

A.j.plus.l 
Epsilon. j 
Epsilon. j.plus.l 
J 


in  Real. Vector; 
in  out  Real.Vector; 
in  out  Real.Vector; 
:  in  integer; 

:  in  integer; 

:  in  out  integer; 

:  in  out  integer; 

:  in  integer)  is 


Pi 

Start .File 

End_File 

Counter 

Window.Boundary .Value 


constant 

integer 

integer 

integer 

constant  lloat 


=  3.141592654; 

=  Data’lirst; 

=  Data 'last; 

=  0; 

=  1.0/sqrt(2.0) ; 
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begin 


- ♦♦♦♦♦CREATE  WIIDOW***** 

—  Window  Defined  Page  11 


—  left  rising  edge 
if  Epsilon. j  *  0  then 

Window. j (A_j)  :=  Window.Boundary .Value; 

else 

for  x  in  (A.j  -  EPSILOI.j) . . (A.j  +  EPSILOI.j)  loop 
Window  j(x)  :=  sin((Pi/(4.0*float(EPSIL0I_j)))  * 

(float (x- C A.j -EPSILOI.j )))) ; 

end  loop; 
end  if; 


—  center 

for  x  in  (A.j  +  EPSILOI.j  +  1) (A.j.plus.l  -  EPSILOI.j .plus. 1  -  1) 
loop 

Window. j(x)  :=  1.0; 
end  loop; 


—  right  falling  edge 
if  EPSILOI.j .plus. 1  *  0  then 

Window. j (A.j.plus.l)  :=  Window.Boundary .Value ; 

else 

for  x  in  (A.j.plus.l  -  EPSILOI.j .plus.l) . . 

(A.j.plus.l  +  EPSILOI.j .plus.l)  loop 
Window. j (x)  :=  cos((Pi/(4.0*float(EPSIL0I_j_plus_l)))  * 

(float (x- (A.j .plus.l  -EPSILOI.j.plus.l)))) ; 

end  loop; 
end  if; 


- ♦♦♦♦♦Multiply  By  Windows**** - 

—  Left  Side 
if  j  *  0  then 

Folded.Data(A.j)  :=  0.0; 
for  x  in  (A.j  +  l)..(A_j  +  EPSILOI.j)  loop 
Folded.Data(x)  :=  Data(x)  *  Window. j(x)  - 

Data(End_File  -  x)  *  Window. j(2*A_j  -  x); 

end  loop; 
else 

for  x  in  (A.j).. (A.j  +  EPSILOI.j)  loop 

Folded.Data(x)  :=  Data(x)  *  Window. j(x)  - 

Data(2*A_j  -  x)  *  Window. j (2*A_j  -  x); 

end  loop; 
end  if; 
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—  Center 

lor  x  in  (A.j  +  EPSILOI.j  +  1) . . U_j_plua_l  -  EPSILOI.j .plus. 1  -  1) 
loop 

Folded.Data(x)  :=  Data(x)  *  Window.j(x); 
end  loop; 

—  Right  Edge 

il  A.j.plus.l  =  Data’laet  then  — (j  =  partitions' last- 1) 

lor  x  in  (A.j.plus.l  -  Epsilon. j.plus.:* )..  (A_j_plus_l  -  1)  loop 
Folded.Data(x)  :=  Data(x)  *  Window.j(x)  + 

Data(A_j_plus_l  -  x)  * 

Window. j(2*A_j_plus_l  -  x); 

end  loop; 

Folded_Data(A_j_plus_l)  :*  2.0*Data(A_j_plus_l)  * 

Window. j (A.j.plus.l) ; 

else 

lor  x  in  (A.j.plus.l  -  EPSILOI.j .plus. 1) .. (A.j.plus.l)  loop 
Folded.Data(x)  :=  Data(x)  *  Window. j(x)  + 

Data(2*A_j_plus_l  -  x)  *  Window. j (2*A_j_plus_l  -  x); 

end  loop; 
end  il; 


end  Multiply .By .Window; 
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—  Compute.Coellicients 

—  This  procedure  takes  in  data  Iron  0  -  I  or  A.j  -  A_j+1  (1+1  points)  and 

—  returns  coellicients  stored  Iron  0  -  1-1  or  A.j  -  (A_j+1  -  1)  (I  points) 

—  Ith  point  is  used  in  calculations  but  is  not  a  valid  out  point 

—  It  will  be  used  to  store  the  zero'th  coellicient  in  the  next  data  segment 


procedure  Compute.Coellicients  (Data.Segment  :  in  out  Real.Vector)  is 

—  Data.Segment  is  Irom  A.j  . .  A.j.plus.l 
Two.I  :  integer  :=  2* (Data.Segment 'last  -  Data.Segment ’lirst) ; 

End.FFT.Data  :  integer  :=  Data.Segment 'lirst  +  Two.I  -  1; 

FFT.Data  :  Complex.Vector (Data.Segment 'lirst. .End.FFT.Data) ; 


begin 


—  (1)  Even  Extend  Data_Segment(See  Vector.Package-Complex.Package) 
FFT.Data  :=  Complex_01(Even_Extend(Data_Segment)) ; 

—  (2)  Perform  Inverse  FFT  (See  FFT.Pack) 
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FFT  (FFT_Data,  True);  —  includes  1/1  scaling 

—  (3. a)  Assign  zero  coefficient  to  first  point  in  Data.Segment 
Data.Segment (Data.Segment ’ first)  :=  FFT_Data(FFT_Data'first) .Real 

*  sqrt  (float  (Two  JO) ; 

—  (3.b)  Assign  coefficients  1  to  1-1  back  to  Data.Segaent 
for  k  in  (Data.Segment’ first  +  1) . . (Data.Segaent' last- 1)  loop 

Data_Segment(k)  :  =  Data.Segaent(k-l)  + 

(2.0  *  sqrt (float (Two JO)  *  FFT.Data(k) .Real) ; 

end  loop; 


end  Compute.Coefficients; 
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—  Reconstruction  FFT 


procedure  Reconstruction_FFT(Input_Data  :  in  Real.Vector; 

Output.Data  out  Real.Vector  )  is 


Pi 

X 

Scale.Factor 

Complex.Factor 

Two.I 

Two.I.Float 

End.FFT.Data 

FFT.Data 


constant 

float 

float ; 

Complex; 

integer 

float 

integer 


3. 141S926B4; 

0.0;  —  float  counter 


2  *  Input .Data’ length; 
float (Two.I) ; 

Input .Data' first  +  Two.*  - 


1; 


Complex.Vector (Input .Data’ first . . End_FFT_Data) ; 


begin 


—  (1)  Odd  Extend  Input.Data  (See  Vector.Package,  Complex.Package) 
FFT.Data  :=  Complex.Of (Odd.Extend(Input.Data) ) ; 


—  (2)  Perform  Inverse  FFT  (See  FFT.Pack) 
FFT(FFT_Data,  True); 

X  :=  0.0; 

Scale.Factor  :=  Bqrt (Two.I.Float ) ; 

—  Range  for  the  scaling  must  be  from  0  to  1-1 
for  1  in  FFT.Data ’range  loop 

Complex.Factor . real  :=  Cos(Pi*X/Two_I_Float) ; 


C-6 


Complex.Factor. imag  :=  Sin(PieX/Two_*_Float) ; 
FFT.Data(l)  :=  Complex.Factor  *  FFT_Data(l); 

X  :=  X  +  1.0;  —  float  counter 

and  loop; 

—  scale  and  assign  imaginary  part  to  Output.Data  1 . .  I 
for  1  in  Output _Data’ range  loop 

Outpnt_Data(l)  ; =  Scale.Factor  *  FFT_Data(l) . imag; 
end  loop; 

end  Reconstruct ion.FFT; 


Divide_By_Windows 


procedure  D ivide_By _H indows  (  Temp.Data 

Data.Segment 
Window .Right 

A-j 

A_j_plus_l 
Epsilon.j 
Epsilon. j.plus.i 
J 


:  in  Real.Vector; 

:  out  Real.Vector; 
:  in  Real.Vector; 
:  in  integer; 

:  in  integer; 

:  in  integer; 

:  in  integer; 

:  in  integer)  is 


Counter 
End_File 
Window .Left 


:  integer  :  =  0; 

:  integer  :=  Temp .Data ’last; 

:  Real.Vector (Vindow.Right 'range) ; 


begin 

Window.Left  :=  Reverse_Assignment(Window_Right) ; 

—  Window.left  is  nos  the  opposite  slope  of  Vindow.Right 

—  and  can  be  used  for  the  j-1  windows  and  Vindow.Right 

—  will  be  used  for  the  jth  windows. 


—  Calculate  Data.Segment  (A.j  -  e.j  +1..  A.j  -  1) 
if  j  =  0  then 

for  i  in  (A.j  -  Epsilon.j  +  1) . . (A.j  -  1)  loop 

Data_Segment(i)  :=  Vindow.Left(i)  *  Temp.Data(End_File  +  i)  - 
Vindow_Left(2*A_j  -  i)  * 

Temp_Data(2*A_j  -  i); 

end  loop; 

else 

for  i  in  (A.j  -  Epsilon.j  +  1) . . (A.j  -  1)  loop 

Data_Segment(i)  :=  Vindow.Left(i)  *  Temp.Data(i)  - 
Vindow_Left(2*A_j  -  i)  * 
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and  loop; 
and  if; 


Te*p._Data(2*A_j  -  i); 


—  Calculate  Data.Segnent  (A.j) 
if  j  =  0  then 

Data.Segnent (A.j)  :=  Teap_Data(End_File)  /  (2.0  *  Vindov.Left(A.j)) ; 

else 

Data.Segnent (A.j )  :=  Tanp_Data(A_j)  /  (2.0  *  Hindov.Left(A.j)) ; 
and  if; 


—  Calculate  Data_Segnent  (A.j  +  1  ..  A_j  +  a_j  -  1) 
if  j  =  0  than 

for  i  in  (A_j  +  1) . . (A_j  +  Epsilon_j  -  1)  loop 
Data.Segnent (i)  :=  Window_Right(2*A_j  -  i)  * 
Tenp_Data(End_File  -  i)  + 
Vindos.Right ( i)  *  Teap_Data(i) ; 

end  loop; 

else 

for  i  in  (A_j  +  1) . . (A_j  +  Epsilon.j  -  1)  loop 
Data_Segnent(i)  :=  Hindoa_Right(2*A_j  -  i)  * 
Temp_Data(2*A_j  -  i)  + 
Uindou.Right(i)  *  Tenp_Data(i) ; 

end  loop; 
end  if; 


—  Calculate  Data_Segment(A_j  +  e_j  . .  A_j-1  -  e_j-l) 

for  i  in  (A_j  +  Epsilon. j) . . (A_j_plus_l  -  Epsilon.j.plus.l)  loop 
Data_Segnent(i)  :=  Tenp.Data(i) ; 
end  loop; 


end  Divide.By .Windows ; 


(Recon.Data 

in 

Real.Vector; 

Actual.Data 

in 

Real. Vector; 

L2_Error 

in 

out 

float; 

RNS.Error 

in 

out 

float; 

L.IlF.Error 

in 

out 

float ; 

lorn.RIV 

in 

out 

float)  is 

error, 
sunl, 
sun2, 
sum3 , 
count 


:  float  :=  0.0; 
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Recon.RIV  :  float  : 3  0.0; 

Actual.RlV  :  float  :=  0.0; 


begin 


for  1  in  Actual_Data' range  loop 

if  abs(Actual.Data(l)  -  Recon.Data(l))  >  orror  then 
error  :=  abs(Actual_Data(l)  -  R«con_Data(l)) ; 
and  if; 

swl  :=  Actual_Data(l)*Actual_Data(l)  +  sunl; 
sua2  :=  Recon_Data(l)*Recon_Data(l)  *  som2; 
sua3  : =  (Actual_Data(l)-Recon_Data(l))* 

(Actual_Data(l)-Rocon_Data(l)>  +  son3 ; 
count  :=  count  +  1.0; 
if  1  =  1  then 
sunl  :=  0.0; 
sun2  :  =  0.0; 

8un3  :=  0.0; 
error : =  0.0; 
count  :*  0.0; 
end  if; 
end  loop; 


RMS.Error 

L2_Error 

L_IIF_Error 

Recon_RIV 

Actual.RlV 

Iorm_RIV 


sqrt (sun3) /count ; 

sqrt (abs(sunl-sum2)/8unl) ; 

error; 

ln(sun2)/ln(2.0) ; 
ln(suml)/ln(2.0); 

Recon_RIV/Actual_RIV  -  1.00000000; 


end  Error _0ut; 
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Do.Work 


procedure  Do.Work  (  Data  in  out  Real. Vector; 

Partitions  :  in  out  Partition.Vector; 

Alphafile.String  :  in  String; 
Outfile.String  :  in  String)  is 


Define  Data  Vectors  and  Their  Ranges. 

FFT.Data  Real.Vector (Data 'range) ; 

Transformed.Data  , 

Temp.Data  Real_Vector((Data’first-l) . . (Data’ last)) ; 

Output.Data  Real.Vector 

((Data’f irst-Partitions(O) .Overlap-1) . . (Data’last+Partitions(0) .Overlap)) ; 
Window  Real.Vector 

((O-Partitions(O) .Overlap) .. (Data 'last  +  Partitions(O) .Overlap)); 
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Start .Data,  End_Data, 

— 

A.j  and  A_j+1 

Start .Window,  End.Window , 

A.j  -  e.j  and  A_j+1  + 

Last .Output .Point , 

— 

A_j+1  -  e_j+l 

Last.Window.Point , 

— 

*-j  +  «-j 

Counter 

integer  : =  0; 

Two.I 

float  :=  0.0; 

Pi 

constant  :=  3.141692664; 

Temp, 

J.Count 

integer  : =  0; 

End.File 

integer  : =  Data’ last; 

Start .File 

integer  :=  Data ’first; 

Alpha! ile, 

Out! ile 

: 

Text.Io . File.Type ; 

L2_Error  : 

float  := 

0.0; 

L.IIF.Error 

float  := 

0.0; 

RMS .Error 

float  : = 

0.0; 

■orm.RIV 

float  : = 

0.0; 

begin 

Text.Io . Create  (Alpha! ile,  Text.Io. Out .File,  Alpha! ile.String) ; 
Text.Io . Create  (Out! ile,  Text.Io. Out .File,  Outfile.String  ); 
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J .Count  : =  -1; 

!or  j  in  0. . (Partitions 'last-1)  loop 
J.Count  : =  J.Count  +  1; 

begin  —  each  segment  is  transformed  within  this  loop 

if  J.Count  =  0  then 

Temp  : =  1;  —  To  set  A(0)  =  0. 
else 

Temp  : =  0; 
end  if; 


- Assign  j  Dependent  Values - 

Start.Data  :=  Partitions(j) .Boundary  -  Temp; 

End.Data  :=  Partitions (j+1) .Boundary ; 

Start.Vindow  :=  Start.Data  -  Partitions(j) .Overlap; 

End. Window  :=  End.Data  +  Part it ions (j+1) .Overlap; 

Two.!  : =  float (2  *  (End.Data-Start.Data)) ; 


Text.Io. put ("Working. . . .  Processing  Data  Segment"); 
Integer_Io.put(j+l) ; 

Text.Io . lew.Line ; 

Mult iply .By _W indow (Data (Data ’ range ) , 

Transformed_Data(Start_Data. .End.Data) , 
Window (Start .Window. . End. Window) ,  Start .data. 


C-10 


End.Data,  Partitions (j) .Overlap, 

Part it ions ( j +1 ) .Over lap,  J .Count) ; 

—  Stops  2  and  3  of  Cosfficisnt  Evaluation  (page  14) 
—  Window  added  as  passed  parameter  only  so 

—  it  does  not  have  to  be  recalculated 

—  in  Divide.By.Vindows 

—  Step  4  of  Coefficient  Evaluation  (page  15) 

Counter  :=  0; 

for  1  in  Start .Data. .End.Data  loop 

Transf ormed.Data(l)  :=  Transf oraed.Data(l)  * 

sin((Pi  *  float (Counter))  /  (Two_I)); 
Counter  :  =  Counter  +1;  —  Counter  goes  from  0  to  1 

end  loop; 


Coopute.Coeff icients  (Transforaed.Dat a (Start .Data. .End.Data)); 

—  This  routine  performs 

—  (1)  an  even  extension  of  the  data, 

—  (2)  an  inverse  FFT, 

—  (3)  assigns  coefficients  0  -  1-1  to  Transf ormed.Data 
—  Steps  S,  6,  and  7  of  Coefficient  Evaluation  (page  16) 

- Print  Alphas  for  each  segment  to  file - 

for  x  in  Start .Data. . ( End.Data- 1)  loop 

f loat.io .put (Alphafile , Transf ormed_Data(x) ) ; 
text_io.new_line(Alphaf ile) ; 
end  loop; 


Reconstruction_FFT(Transformed_Data( (Start .Data) . .  (End_Data-l)) , 
FFT_Data( (Start _Data+l) . .End_Data)) ; 

—  This  routine  performs 

—  (1)  an  odd  extension  of  the  data, 

—  (2)  an  inverse  FFT, 

—  (3)  assigns  the  scaled  imaginary  result  to 

the  Output.Data  values  1  -  I  or  A_j  +  1  to  A_j+1 
—  This  indicates  that  data  values  0-1 

—  create  alphas  from  0  -  1-1  which  in  turn  are  used 

—  to  reconstruct  data  values  from  1  -  I 

—  Steps  1  and  2  of  Reconstruction  (page  15) 


end; - block 

end  loop;  —  —  j  loop 
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J.Count  : =  -1; 

for  j  in  0. . (Partitions ’last-1)  loop 
J.Count  :=  J.Count  +  1; 
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begin  —  etch  swgnwnt  is  tranaforaed  within  this  loop 

if  J.Count  =  o  then 
Tsap  :  =  1; 
slss 

Tsap  :  =  0; 
end  if; 

Assign  j  Ospsndsnt  Values - 

Start.Data  :  =  Partitions(j)  .Boundary  -  Top; 

End. Data  :*  Part it ions (j+1) .Boundary; 

Start.Window  :*  Start.Data  -  PartitionsCj) .Overlap; 

End_ Window  :*  End.Data  +  Part it ions (j+1) .Overlap; 

Two-1  ••=  float  (2  *  (End_Data-Start_Data) ) ; 

Last _0ut put .Point  :=  End.Data  -  Partitions(j+l) .Overlap; 
Laat.Window.Point  :=  Start.Data  +  Partitions(j) .Overlap; 


Divide. By. Windows (FFT.Dat a (FFT.Data ' range) , 

Output _Data( (Start_Window+l) . . Las t.Output .Point) , 
Window (Start .Window. . Last .Window .Point) , 

Start .Data,  End_Data,Partitions(j) -Overlap, 

Part it ions (j+1). Over lap,  J.Count); 

—  Move  the  (0. .  (-Epsilon(0)+1))  Data  Back  To  The 
—  (End.File. . (End_File-Epeilon(0)+l))  Range, 
if  j .count  a  0  then 

for  k  in  0. . (Partitions(O) .Overlap-1)  loop 

Output _Data(End_File  -  k)  :=  Output .Data(-k) ; 
end  loop; 
end  if; 


end;  —  —  block 
end  loop;  —  —  j  loop 


- CALCULATE  OUTPUT  ERRORS  - 

Error.Out  (  Output .Data(Data’range) , 

Data (Data’ range) , 

L2_Error ,  RMS.Brror,  L.IIF.Error,  Bom.RIV) ; 


-  PRIBT  OUT  TO  FILES  - 

Text.io . Put .Line ("print  output  to  file"); 
Text.Io . lew.Line ; 

Text.Io. put ("L2  Error  =  "); 
Float_Io.put(L2_Error) ; 

Text.Io . lew.Line ; 
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Text.Io.  put  (**L(  Infinity)  Error  ■  "); 

Float _Io .put (L_IIF_Error) ; 

Text.Io . les.line ; 

Text.Io. pat (“RMS  Error  »  "); 

Float_Io . put (RKS_ Error ) ; 

Text.Io . lev.Line ; 

Text.Io. put("Iormalized  RIV  =  H); 

Float _Io .put (Ior»_RIV) ; 

Taxt.Xo . lev.Line ; 

for  x  in  Data ’range  loop 

Float _Io . put (Outf ila , ( ( Output _Dat a (x) ) ) ) ; 
Text.Io . lev.Line (Outf ila) ; 

Text.Io. lev_Line(Outf ila) ; 
and  loop; 

Taxt.Xo . Cloaa(llphaf ila) ; 
Taxt.Io.Cloaa(Outfila) ; 


and  Do.Vork; 


begin  — sain 

Gat.Fila.Data  (Alpha_Iaae,  ilpha.Iaaa.Langth , 

Outf ile.Iaae,  Outf ile.Iaae .Length, 
Big_Daddy_Partition_Vector ,  Actual.Partitions , 
Big_Daddy_Data_Vector,  Actual.Data.Langth) ; 

—  By  sanding  the  correct  size  arrays  into  Do.Vork,  the  array 

—  attributes  can  ba  us ad  to  determine  the  size  rather  than 

—  passing  another  pax  ana tar 

Do.Vork  (Big_Daddy_Data_Vector ( 1 . . Actual_Data.Length) , 

Big_Daddy_Partition_Vector  (0. .Actual.Partitions), 
Alpha_laae(l . . Alpha.Iane.Length) , 

Outfila.laaad.  .Outfila.Iaaa.Langth)) ; 


and;  — aain 
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Appendix  D.  Ada  Source  Code  For  Linear  Transforms 


The  following  code  was  written  to  perform  the  linear  transforms  for  the  homomorphic  and 
non-homomorphic  systems  studied  in  Chapter  IV. 


-  Linear  Trans fora  - 

procedure  Linear .Transform (  Inpat .Data  :  in  Eeal.Vector; 

Output .Data  :  out  Eeal.Vector; 
Inverse  :  in  out  Integer; 

Transform  :  in  out  Integer; 

Min  :  in  out  Eeal.Vector; 

Sign  :  in  out  Eeal.Vector; 

J  :  in  Integer  )  is 

x  :  integer; 
teap  :  float; 

mult  :  float  :  =  10.0000; 
add  :  float  :=  1.0000; 

begin 

if  Inverse  =  0  then 

— Find  minimus  value  for  each  uindov 
for  1  in  Input .Data’ range  loop 
if  Input_Data(l)  <  min(J)  then 
Nin(j)  :=  Input _Dat a(l) ; 
end  if; 
end  loop; 

if  Transform  =  1  then 

for  1  in  Input .Data’ range  loop 
Output.Data(l)  :  =  Input.Data(l) ; 
end  loop; 
end  if; 

if  Transform  =  2  then 

for  1  in  Input .Data’ range  loop 
if  Input _Data(l)  <0.0  then 
sign(l)  :=  -1.0; 

else 

sign(l)  :=  1.0; 
end  if; 

Output.Data(l)  :=  sqrtfabs (Input _Data(l))) ; 
end  loop; 
end  if; 

if  Transform  =  3  then 

for  1  in  Input .data’ range  loop 
if  Input_Data(l)  <0.0  then 
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sign(l)  :=  -1.0; 

•Isa 

•ign(l)  :=  1.0; 

•nd  if; 

Output.Data(l)  :=  ln(abs ( Input _Data(l))) ; 

•nd  loop; 

•nd  if; 

if  Trans fora  =  4  than 

for  1  in  Input.data’rang*  loop 

Ontput_Data(l)  :»  anlt  *  Input _Data(l)  +  add; 

•nd  loop; 
end  if; 

if  Transform  =  5  than 

for  1  in  Input .Data’ rang*  loop 
Output  JData(l)  : =  In  (  Input _Data(l)  ♦ 

sqrt(Input_Data(l)*Input_Data(l)  +  1.0)  ); 

•nd  loop; 

•nd  if; 

•nd  if;  —  End  Forward  Transfora  Loop 


if  Inverse  =  1  then 

if  Transform  =  1  than 

for  1  in  Input .Data’ rang*  loop 

Output _Data(l)  :=  Input _Data(l) ; 

•nd  loop; 

end  if; 

if  Transform  =  2  than 

for  1  in  Input_Oata'ranga  loop 

taap  :=  Input _Data(l)*Input_Data(l) ; 
if  sign(l)  =  -1.0  than 

Output_Data(l)  :=  -taap; 

•lse 

Output JData(l)  :=  taap; 

•nd  if; 

•nd  loop; 

•nd  if; 

if  Transform  =  3  than 

for  1  in  Input .Data’ range  loop 
taap  : =  exp(Input_Data(l)) ; 
if  8ign(l)  =  -1.0  then 
Output _Data(l)  : =  -temp; 

•lse 

0utput_Data(l)  :=  taap; 

•nd  if; 
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end  loop; 
end  if; 

if  Tranefora  =  4  then 

for  1  in  Input _Data’ range  loop 

Output _Data(l)  : =  ( Input _Data(l)  -  add)/ault; 
end  loop; 
end  if; 

if  Tranafora  =  8  then 

for  1  in  Input _data» range  loop 
Output _Data(l)  :=  0.80  *  (  exp(Input_Data(l))  - 

exp(-Input_Data(l))  ); 

end  loop; 
end  if; 

end  if;  —  End  Inverse  Transfora  Loop 


end  Linear.Transfora; 
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Appendix  E.  Ada  Source  Code  For  Compression  and  Error  Calculations 


Percent  Compression 

The  following  source  code  was  used  to  perform  the  compression  of  coefficients  in  the  experi¬ 
ments  of  Chapters  III,  IV,  and  V.  It  is  in  a  generic  form,  because  the  ranges  for  the  loops  of  the 
Compression  package  are  determined  by  what  precedes  it  (complex  MW,  sinusoidal  MW,  or  one  of 
the  systems  presented  in  Chapter  IV).  The  data  that  is  changed  (in  this  version  “Output-Data”) 
also  depends  upon  what  precedes  it  (for  example,  the  complex  MW  uses  “FFT_Data”  at  this  point 
of  the  program. 


- Compress  the  Homomorphic  Coefficients  - 

Io_Zeros_R  :=  float (Percent ..Compression)  * 
float (data ' last) ; 
lo.Zeros  :=  integer (*o_Zeros_R) ; 

Text_Io . Mem.line ; 

Text.Io. put ("lo.Zeros  =  "); 

Integer_Io.put(*o_Zeros) ; 

Text.Io . lew. line ; 

—  Find  max  output  to  establish  initial  threshold  and  to  set  the 

—  level  for  adjustment, 
for  x  in  (  )  . .  (  )  loop 

if  Output_Data(x)  >  Max  then 
Max  :=  Output _Data(x) ; 
end  if; 
end  loop; 

Cep.Threshold  : =  Percent .Compression  *  Max; 

Adjust  :=  0.10  *  Max; 


loop  —  Start  Infinite  Loop 

for  x  in  (  )  . .  (  )  loop 

—  Compression  Output  for  each  segment  to  file- 

if  abs (Output _Data(x))  <  Cep.Threshold  then 
Temp(x)  :=  0.0; 

Compression.Count  : =  Compression.Count  +  1; 
else  Temp(x)  :=  1.0; 
end  if; 
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«nd  loop;  —  the  for-loop  ins id*  of  the  infinite  loop 


if  (Compression_Count  <  Io_*eros)  then 
if  TooLov  >  0  than 

Cep.Threshold  : =  Cep_Threshold  +  Adjust; 
end  if; 

if  TooHigh  >  0  thsn 

Adjust  :=  Adjust/2.0; 

Cep_ Threshold  :=  Cep_Threshold  +  Adjust; 

•nd  if; 

TooLov  : =  1 ; 

TooHigh  :=  0; 
end  if; 

if  (Compression_Count  >  Io_Zeros)  then 
if  TooHigh  >  0  then 

Cep.Threshold  :=  Cep_Threshold  -  Adjust; 
end  if; 

if  TooLov  >  0  then 

Adjust  :=  Adjust/2.0; 

Cep.Threshold  :=  Cep_Threshold  -  Adjust; 
end  if; 

TooHigh  :=  1; 

TooLov  : *  0 ; 
end  if; 

if  (Compress ion_Count  =  lo.Zeros)  then 
exit; 
end  if; 

Compression.Count  :=  0; 

Cep_Threshold  : =  Cep_Threshold; 

Adjust  :*  Adjust; 

end  loop;  — end  of  infinite  loop 

for  x  in  (  )  . .  (  )  loop 

—  Set  compressed  Output  values  to  zero 
Output  J)ata(x)  :=  Output _Data(x)  *  Temp(x); 

—  Print  output  values  to  the  alpha  file  (coefficients) . 
Float _Io. Put  (Alphafile,  Output _Data(x)) ; 

Text_Io . Hev_Line (Alphafile) ; 
end  loop; 


Error  Calculations  Code 

The  error  calculation  code  below  is  for  all  real  valued  outputs  and  inputs.  The  complex  MW 
produces  complex  outputs.  The  user  must  ensure  that  the  correct  type  of  data  is  being  sent  to  this 


package.  It  was  assumed  that  all  inputs  would  be  red  and  therefore  all  outputs  from  the  CMW 
would  also  be  real.  In  this  case  only  the  real  part  of  “Recon_Data”  is  used.  Also  the  RIV  criteria 
that  is  calculated  is  actually  (1  -  RIV),  because  the  values  are  often  so  close  to  one  that  it  is  difficult 
to  distinguish  them. 


- ERROR.OUT - 

procedure  Error.Out  (Recon.Data  :  in  Real_Vector; 

Actual.Data  :  in  Real. Vector; 
L2_Error  :  in  out  float; 

RMS.Error  :  in  out  float; 

L.IIF.Error  :  in  out  float; 

lorm.RIV  :  in  out  float)  is 


error , 

sunl , 

8um2 , 
sum3, 
count 

:  float 

:=  0.0; 

Recon.RIV 

:  float 

:=  0.0; 

Actual.RIV 

:  float 

:=  0.0; 

begin 

for  1  in  Actual.Data’ range  loop 

if  abs(Actual_Data(l)  -  Recon.Data(l))  >  error  then 
error  :=  abs(Actual_Data(l)  -  Recon.Data(l)) ; 
end  if; 

suml  :=  Actual_Data(l)*Actual_Data(l)  +  suml; 
sum2  :=  Recon_Data(l)*Recon_Data(l)  +  sum2; 
sun3  :=  (Actual.Data(l)-Recon.Data(l))* 

(Actual.Data(l)-Recon.Data(l))  +  sum3; 
count  : =  count  +  1.0; 
end  loop; 

RMS.Error  :  =  sqrt(sum3) /count; 

L2_Error  : =  sqrt(abs(sual-sum2)/sual) ; 

L.IIF.Error  :=  error; 

Recon_RIV  : =  ln(sum2)/ln(2.0) ; 

Actual.RIV  :=  ln(sual)/ln(2.0) ; 

lorm.RIV  :=  Recon.RIV/Actual.RIV  -  1.00000000; 

end  error.out 
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